#ifndef OBSERVABLES_H
#define OBSERVABLES_H
#include "utility.h"
#include <stdio.h>
#include <iostream>
#include "mt2_bisect.h"
#include "TRandom1.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include "TStopwatch.h"
#include "TF1.h"
#include "TUnuran.h"
#include "TUnuranContDist.h"

#define PI 3.1415926535897932384626433832795

using namespace std;
using namespace mt2_bisect;

struct CompositeParticle {
	TLorentzVector Particle;
	vector<TLorentzVector> Components;
};

//****************************************************************************
extern "C" long double beta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double beta;

	if (vectorstoadd.size() == 1) {
		beta = vec_part[vectorstoadd[0]].Beta();
	} else {
		cout << " beta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return beta;
}

//****************************************************************************
extern "C" long double betaofn_sign(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	res = sum.Pz() / sum.E();

	return res;
}

//****************************************************************************
extern "C" long double betaofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	res = sum.Beta();

	return res;
}

//****************************************************************************
extern "C" long double lorentzgamma(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double gamma;

	if (vectorstoadd.size() == 1) {
		gamma = vec_part[vectorstoadd[0]].Gamma();
	} else {
		cout << " gamma can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return gamma;
}

//****************************************************************************
extern "C" long double lorentzgammaofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	res = sum.Gamma();

	return res;
}

//****************************************************************************
extern "C" long double phi(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double phi;

	if (vectorstoadd.size() == 1) {
		phi = vec_part[vectorstoadd[0]].Phi();
	} else {
		cout << " phi can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return phi;
}

//****************************************************************************
extern "C" long double phiofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double phiofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	phiofn = sum.Phi();

	return phiofn;
}

//****************************************************************************
extern "C" long double thrust(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, thisT, norm, theta, divisions;
	long double px, py;

	res = 0;
	theta = 0;
	divisions = 2000.;

	for (int step = 0; theta < 2 * PI; step++) {
		norm = 0;
		thisT = 0;

		theta = theta + 2. * PI / divisions;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			px = vec_part[vectorstoadd[vect]].Px();
			py = vec_part[vectorstoadd[vect]].Py();

			thisT = thisT + fabs((px * cos(theta) + py * sin(theta)));

			norm = norm + sqrt(px * px + py * py);

			//cout << "theta is currently " << theta << " and T=" << thisT << " so far"<< endl;
		}
		//cout << "for theta " << theta << " the thrust is " << thisT / norm
		//<< endl;
		if (thisT / norm > res) {
			res = thisT / norm;
			//cout << "*** now the thrust is " << res << endl;
		}

	}

	return res;
}

//****************************************************************************
extern "C" long double pTaverageofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

//	printvector(vectorstoadd);

	long double result;
	long double HT = 0;
	int vect = 0;

	long double numberofparticles = vectorstoadd.size();

	while (vect < vectorstoadd.size()) {
		HT = HT + fabs(vec_part[vectorstoadd[vect]].Pt());
		vect++;
	}
	result = HT / numberofparticles;
//	cout << "the number of paticles is " <<  vectorstoadd.size() << endl;
//	cout << "the observable is " <<  result << endl;
	return result;
}

//****************************************************************************
extern "C" long double HT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double HT = 0;
	int vect = 0;

	while (vect < vectorstoadd.size()) {
		HT = HT + fabs(vec_part[vectorstoadd[vect]].Pt());
		vect++;
	}

	return HT;
}

//****************************************************************************
extern "C" long double muTH(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	float beta = 1;

	long double muQ = 81.0 * 81.0;
	int vect = 0;

	while (vect < vectorstoadd.size()) {

		muQ = muQ
				+ 0.5 * fabs(vec_part[vectorstoadd[vect]].Pt())
						* fabs(vec_part[vectorstoadd[vect]].Pt());
		vect++;
	}

	return beta * sqrt(muQ);
}

extern "C" long double Meffofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// this is a property of a group of particles in relation with the mET
	// bbar.exe Meffofn 55,-55,211234 0
	// then one can select a particular set of jets ordered by pT

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double meff, res;

	int vect = 0;
	TLorentzVector sum, visible;
	vector<int> neutrinos;
	neutrinos.push_back(121416);

	// scan all particles in the event and find the mET from Delphes
	for (int vect = 0; vect < part.size(); vect++) {

//		if ((part[vectorstoadd[vect]]->Status > 0)
//			&& (part[vectorstoadd[vect]]->Status != 2)) {

//		cout << " Particle " << vect << " " << vec_part[vect].Pt() << " PID: " << part[vect]->PID <<  endl;

		if (isinclass(neutrinos, part[vect]->PID) == true) {

			// it is invisible, thus sum with mET
			sum = sum + vec_part[vect];
		}

//		}

	}

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		visible = visible + vec_part[vect];
	}

	meff = sum.Pt() + visible.Pt();

	//cout<<pTofn<<endl;
	return meff;
}

extern "C" long double f(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// this is a property of a group of particles in relation with the mET.
	// the particles that make the mET are hard-coded here as 12,14,16,-12,-14,-16 and 121416
	// Usage:  bbar.exe f 55,-55,5,-5 0
	// then one can select a particular set of jets ordered by pT in the single particles filters

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double meff, f, HT;

	HT = 0;
	int vect = 0;
	TLorentzVector sum;
	vector<int> neutrinos;
	neutrinos.push_back(121416);
	neutrinos.push_back(1000022);
	neutrinos.push_back(12);
	neutrinos.push_back(-12);
	neutrinos.push_back(14);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	neutrinos.push_back(16);

	// scan all particles in the event and find the mET from Delphes
	for (int vect = 0; vect < part.size(); vect++) {

//		if ((part[vectorstoadd[vect]]->Status > 0)
//			&& (part[vectorstoadd[vect]]->Status != 2)) {

//		cout << " Particle " << vect << " " << vec_part[vect].Pt() << " PID: " << part[vect]->PID <<  endl;

		if (isinclass(neutrinos, part[vect]->PID) == true) {

			// it is invisible, thus sum with mET
			sum = sum + vec_part[vect];
		}

//		}

	}

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		HT = HT + fabs(vec_part[vectorstoadd[vect]].Pt());
	}

	meff = sum.Pt() + HT;

	f = sum.Pt() / meff;
	//cout<<pTofn<<endl;
	return f;
}

extern "C" long double mT_ofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// this finds the mET in the delphes event and computes the transverse mass of the
	// visible particles that it has been fed with

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	int vect = 0;
	TLorentzVector mET, visible;
	vector<int> neutrinos;
	neutrinos.push_back(121416);
	neutrinos.push_back(1000022);
	neutrinos.push_back(12);
	neutrinos.push_back(-12);
	neutrinos.push_back(14);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	neutrinos.push_back(16);

	// scan all particles in the event and find the mET from Delphes
	for (int vect = 0; vect < part.size(); vect++) {

//		if ((part[vectorstoadd[vect]]->Status > 0)
//			&& (part[vectorstoadd[vect]]->Status != 2)) {

//		cout << " Particle " << vect << " " << vec_part[vect].Pt() << " PID: " << part[vect]->PID <<  endl;

		if (isinclass(neutrinos, part[vect]->PID) == true) {

			// it is invisible, thus sum with mET
			mET = mET + vec_part[vect];
		}

//		}

	}

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
//		cout<< part[vect]->PID << " pX=" << vec_part[vectorstoadd[vect]].Px() << endl;
		visible = visible + vec_part[vectorstoadd[vect]];
	}

	// ****** Compute mT *******
	// as in the ATLAS TDR pg 706 because I want to reproduce their plot
	float pxV, pyV, pxmiss, pymiss;

	pxmiss = mET.Px();
	pymiss = mET.Py();

	pxV = visible.Px();
	pyV = visible.Py();

	float abs_pTV, abs_pTmiss, pTVdotpTmiss;

	abs_pTV = sqrt(pxV * pxV + pyV * pyV);
	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);
	pTVdotpTmiss = pxmiss * pxV + pymiss * pyV;
	res = sqrt(2 * (abs_pTmiss * abs_pTV - pTVdotpTmiss));

	//cout << "pxV:"<< pxV  << "pyV:"<< pyV  << " pxmiss:" << pxmiss << " pymiss:" << pymiss << " res:"<< res << endl;
	return res;
}

extern "C" long double delta_phi_mETdelphes_single(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// this is a property of a sinlge jet
	// bbar.exe delta_phi_mETdelphes_single 55,-55,211234 1
	// or one can select a particular jets ordered by pT

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double dphi, res;

	int vect = 0;
	TLorentzVector sum, visible;
	vector<int> neutrinos;
	neutrinos.push_back(121416);
	neutrinos.push_back(1000022);
	neutrinos.push_back(12);
	neutrinos.push_back(-12);
	neutrinos.push_back(14);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	neutrinos.push_back(16);

	// scan all particles in the event and find the mET from Delphes
	for (int vect = 0; vect < part.size(); vect++) {

//		if ((part[vectorstoadd[vect]]->Status > 0)
//			&& (part[vectorstoadd[vect]]->Status != 2)) {

//		cout << " Particle " << vect << " " << vec_part[vect].Pt() << " PID: " << part[vect]->PID <<  endl;

		if (isinclass(neutrinos, part[vect]->PID) == true) {

			// it is invisible, thus sum with mET
			sum = sum + vec_part[vect];
		}

//		}

	}

	if (vectorstoadd.size() == 1) {
		TLorentzVector object;
		object = vec_part[vectorstoadd[0]];
		//cout << part[vectorstoadd[0]]->PID << endl;
		//cout << vec_part[vectorstoadd[0]].Pt() << endl;
		dphi = fabs(sum.Vect().DeltaPhi(object.Vect()));
	} else {
		cout
				<< "delta_phi with the mET can be evaluated only for a single particle"
				<< endl;
		exit(1);
	}

	//cout<<pTofn<<endl;
	return dphi;
}

//****************************************************************************
extern "C" long double mET_SM_nosmear(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	long double HT;
	int vect = 0;
	TLorentzVector sum, visible;
	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	//cout << vect << " " << part[vect]->PID << endl;

	//	float length;
	//	length=sizeof(vec_part[8]);
	//	cout<<length<<endl;
	//  int pid;
	//  pid = part[0]->PID;
	//  cout << pid;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		//warning MG showered files contains some status 2 particle, which is not present in pythia output or in fastjet output

		if ((part[vectorstoadd[vect]]->Status > 0)
				&& (part[vectorstoadd[vect]]->Status != 2)) {

			if (isinclass(neutrinos, part[vectorstoadd[vect]]->PID) == true) {

				// it is invisible, thus sum with mET
				sum = sum + vec_part[vectorstoadd[vect]];
			} else if (isinclass(neutrinos, part[vectorstoadd[vect]]->PID)
					== false) {

				// it is visible, thus sum to compute HT which is relevant for the smearing
				HT = HT + fabs(vec_part[vectorstoadd[vect]].Pt());
			}
		}

	}

	pTofn = sum.Pt();

	//cout<<pTofn<<endl;
	return pTofn;
}

//****************************************************************************
extern "C" long double mET_SM_smeared(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn, smear, result, HT;
	Double_t sigma;
	HT = 0;
	TLorentzVector sum, visible;
	vector<int> neutrinos;

	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);

	//	float length;
	//	length=sizeof(vec_part[8]);
	//	cout<<length<<endl;

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		//warning MG showered files contains some status 2 particle, which is not present in pythia output or in fastjet output

		if ((part[vectorstoadd[vect]]->Status > 0)
				&& (part[vectorstoadd[vect]]->Status != 2)) {
			//			puts("found final");

			if (isinclass(neutrinos, part[vectorstoadd[vect]]->PID) == true) {
				// it is invisible, thus sum with mET
				//				puts("found invisible");
				sum = sum + vec_part[vectorstoadd[vect]];
			} else if (isinclass(neutrinos, part[vectorstoadd[vect]]->PID)
					== false) {
				//				puts("not a nu");
				// it is visible, thus sum to compute HT which is relevant for the smearing
				HT = HT + fabs(vec_part[vectorstoadd[vect]].Pt());

			}
		}

	}

	pTofn = sum.Pt();
	//cout << HT << endl;
	sigma = 0.5 * sqrt(HT);
	TRandom *r3 = new TRandom3(0);
	smear = r3->Gaus(0, sigma);
	//	cout<<pTofn<<" "<< HT << " "<<sigma<<" "<<smear<<endl;
	result = pTofn + smear;

	//	cout<<result<<endl;
	return result;
}

//****************************************************************************
extern "C" long double mET(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	pTofn = sum.Pt();
	//cout<<pTofn<<endl;
	return pTofn;
}

//****************************************************************************
extern "C" long double pXofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	pTofn = sum.Px();
	//cout<<pTofn<<endl;
	return pTofn;
}

//****************************************************************************
extern "C" long double pYofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	pTofn = sum.Py();
	//cout<<pTofn<<endl;
	return pTofn;
}

//****************************************************************************
extern "C" long double pZofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	pTofn = sum.Pz();
	//cout<<pTofn<<endl;
	return pTofn;
}

//****************************************************************************
extern "C" long double pTratio(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTa, pTb, ratio;
	int vect = 0;
	TLorentzVector sum;
	if (vectorstoadd.size() == 2) {
		pTa = vec_part[vectorstoadd[0]].Pt();
		pTb = vec_part[vectorstoadd[1]].Pt();

	}

	if (pTa > pTb) {
		ratio = pTb / pTa;
	}

	if (pTa < pTb) {
		ratio = pTa / pTb;
	}

	return ratio;
}

//****************************************************************************
extern "C" long double pTofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	pTofn = sum.Pt();
	/*cout<<pTofn<<endl;*/
	return pTofn;
}

//****************************************************************************
extern "C" long double abs3Vofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double result;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	result = sqrt(
			sum.Px() * sum.Px() + sum.Py() * sum.Py() + sum.Pz() * sum.Pz());
	//cout<<pTofn<<endl;
	return result;
}

//****************************************************************************
extern "C" long double DeltapT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 2) {

		long double pTofn;
		int vect = 0;
		TLorentzVector diff;

		diff = vec_part[vectorstoadd[0]] - vec_part[vectorstoadd[1]];

		//printlv(sum);

		pTofn = fabs(diff.Pt());
		//pTofn = sqrt(diff.Px()*diff.Px() + diff.Py()*diff.Py()  );
		//cout<<pTofn<<endl;
		return pTofn;

	} else {
		cout << "DeltapT can be computed only on a pair" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double deltapT2(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 2) {

		long double pTofn, diffx, diffy, res;
		int vect = 0;
		TLorentzVector diff;

		diffx = vec_part[vectorstoadd[0]].Px() - vec_part[vectorstoadd[1]].Px();
		diffy = vec_part[vectorstoadd[0]].Py() - vec_part[vectorstoadd[1]].Py();
		res = sqrt(diffx * diffx + diffy * diffy);

		return res;

	} else {
		cout << "DeltapT can be computed only on a pair" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double pTminofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTminofn;
	int vect = 0;
	pTminofn = vec_part[vectorstoadd[0]].Pt();

	while (vect < vectorstoadd.size()) {

		long double thispt;
		thispt = vec_part[vectorstoadd[vect]].Pt();

		if (thispt < pTminofn) {

			pTminofn = thispt;

		}

		vect++;
	}

	return pTminofn;
}

//****************************************************************************
extern "C" long double pTmaxofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTmaxofn;
	int vect = 0;
	pTmaxofn = vec_part[vectorstoadd[0]].Pt();

	while (vect < vectorstoadd.size()) {

		long double thispt;
		thispt = vec_part[vectorstoadd[vect]].Pt();

		if (thispt > pTmaxofn) {

			pTmaxofn = thispt;

		}

		vect++;
	}

	return pTmaxofn;
}

//****************************************************************************
extern "C" long double pTminoverminv(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pTminofn, mll, result;

	pTminofn = vec_part[vectorstoadd[0]].Pt();

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		long double thispt;
		thispt = vec_part[vectorstoadd[vect]].Pt();

		if (thispt < pTminofn) {

			pTminofn = thispt;

		}

	}

	TLorentzVector sum;

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		sum = sum + vec_part[vectorstoadd[vect]];
	}
	mll = sqrt(sum * sum);

	result = pTminofn / mll;
	return result;
}

//****************************************************************************
extern "C" long double Eofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double Eofn;
	int vect = 0;
	TLorentzVector sum;
	//puts("involved particles ");
	//printvector(vectorstoadd);
	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	Eofn = sum.E();
	//cout<<pTofn<<endl;
	return Eofn;
}

//****************************************************************************
extern "C" long double ET(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pt;

	if (vectorstoadd.size() == 1) {
		pt = vec_part[vectorstoadd[0]].Et();
	} else {
		cout << " pT can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return pt;
}

//****************************************************************************
extern "C" long double QbranchingNoTag(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	TLorentzVector diff;
	long double pl, M;
	int pid;

	if (vectorstoadd.size() == 1) {

		long double eta, pT, energy;

		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

		if (!(energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5)) {

			pl = vec_part[vectorstoadd[0]].Pz();
			pid = part[vectorstoadd[0]]->PID;

			for (int istate = 0; istate < part.size(); istate++) {
				// loop on the initial states
				if (part[istate]->PID != pid && part[istate]->Status == -1) {
					//check the sign of the rapidity
					if (pl * vec_part[istate].Pz() > 0) {
						diff = vec_part[vectorstoadd[0]] - vec_part[istate];
						M = sqrt(fabs(diff * diff));
						return M;

					}
				}
			}

			return -200;

		} else {
			return -100;
		}

	} else {
		cout << " Q-branching can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double QbranchingTagW(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	TLorentzVector diff;
	long double pl, M;
	int pid;

	if (vectorstoadd.size() == 1) {

		long double eta, pT, energy;

		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

		if ((energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5)) {

			pl = vec_part[vectorstoadd[0]].Pz();
			pid = part[vectorstoadd[0]]->PID;

			for (int istate = 0; istate < part.size(); istate++) {
				// loop on the initial states
				if (part[istate]->PID != pid && part[istate]->Status == -1) {
					//check the sign of the rapidity
					if (pl * vec_part[istate].Pz() > 0) {
						diff = vec_part[vectorstoadd[0]] - vec_part[istate];
						M = sqrt(fabs(diff * diff));
						return M;

					}
				}
			}

			return -200;

		} else {
			return -100;
		}

	} else {
		cout << " Q-branching can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double QbranchingTag(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	TLorentzVector diff;
	long double pl, M;

	if (vectorstoadd.size() == 1) {

		long double eta, pT, energy;

		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

		if (energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5) {

			pl = vec_part[vectorstoadd[0]].Pz();

			for (int istate = 0; istate < part.size(); istate++) {
				// loop on the initial states
				if (part[istate]->Status == -1) {
					//check the sign of the rapidity
					if (pl * vec_part[istate].Pz() > 0) {
						diff = vec_part[vectorstoadd[0]] - vec_part[istate];
						M = sqrt(fabs(diff * diff));
						return M;

					}
				}
			}
		} else {
			return -100;
		}

	} else {
		cout << " Q-branching can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double pTtATLAS(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	TVector3 pT1, pT2, t, sum, pTt;

	long double pTtAbs;

	if (vectorstoadd.size() == 2) {

		pT1 = vec_part[vectorstoadd[0]].Vect();
		pT2 = vec_part[vectorstoadd[1]].Vect();

		pT1.SetZ(0);
		pT2.SetZ(0);

//		cout << "{" << pT1.Px() << ", " << pT1.Py() << ", " << pT1.Pz() << "}" << endl;
		//	cout << "{"<< pT2.Px() << ", " << pT2.Py() << ", " << pT2.Pz() << "}" << endl;

		t = (pT1 - pT2);

		//cout << t.Px() << " " << t.Py() << " " << t.Pz() << endl;

		t.SetMag(1); // keeping theta and phi

		//cout << t.Px() << " " << t.Py() << " " << t.Pz() << endl;

		sum = pT1 + pT2;

//		cout << sum.Px() << " " << sum.Py() << " " << sum.Pz() << endl;

		pTt = t.Cross(sum);

		pTtAbs = pTt.Mag();

		//	cout << pTtAbs << endl;
		//cout << "r37289rwkjehfsjgfsjdhfgiuwrt8e7wryuwiryiwerhiu" << endl;

	} else {
		cout << "pTtATLAS can be evaluate only for 2 fourvectors" << endl;
		exit(1);
	}

	return pTtAbs; //} else {

	//}
}

//****************************************************************************
extern "C" long double pT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pt;

	if (vectorstoadd.size() == 1) {
		pt = vec_part[vectorstoadd[0]].Pt();
	} else {
		cout << " pT can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return pt;
}

//****************************************************************************
extern "C" long double pT3(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pt;
	//cout<<part[vectorstoadd[0]]->PID<<endl;
	if (vectorstoadd.size() == 1) {
		pt = vec_part[vectorstoadd[0]].Pt();
	} else {
		cout << " pT can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	return pt;
}

//****************************************************************************
extern "C" long double betarazor(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double betaR;
	long double Ea, Eb, pZa, pZb;
	if (vectorstoadd.size() == 2) {
		Ea = vec_part[vectorstoadd[0]].E();
		pZa = vec_part[vectorstoadd[0]].Pz();
		Eb = vec_part[vectorstoadd[1]].E();
		pZb = vec_part[vectorstoadd[1]].Pz();
	} else {
		cout << " betarazor can be evaluated only for two particles" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	betaR = (Ea - Eb) / (pZa - pZb);
	return betaR;
}

//****************************************************************************
extern "C" long double MRazorof2(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double MR;
	long double Ea, Eb, pZa, pZb;
	if (vectorstoadd.size() == 2) {
		Ea = vec_part[vectorstoadd[0]].E();
		pZa = vec_part[vectorstoadd[0]].Pz();
		Eb = vec_part[vectorstoadd[1]].E();
		pZb = vec_part[vectorstoadd[1]].Pz();
	} else {
		cout << " betarazor can be evaluated only for two particles" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	MR =
			2
					* sqrt(
							(Ea * pZb - Eb * pZa) * (Ea * pZb - Eb * pZa)
									/ ((pZa - pZb) * (pZa - pZb)
											- (Ea - Eb) * (Ea - Eb)));
	return MR;
}

//****************************************************************************
extern "C" long double MRversusEgeomavg(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double MR, result;
	long double Ea, Eb, pZa, pZb;
	if (vectorstoadd.size() == 2) {
		Ea = vec_part[vectorstoadd[0]].E();
		pZa = vec_part[vectorstoadd[0]].Pz();
		Eb = vec_part[vectorstoadd[1]].E();
		pZb = vec_part[vectorstoadd[1]].Pz();
	} else {
		cout << " betarazor can be evaluated only for two particles" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	MR =
			2
					* sqrt(
							(Ea * pZb - Eb * pZa) * (Ea * pZb - Eb * pZa)
									/ ((pZa - pZb) * (pZa - pZb)
											- (Ea - Eb) * (Ea - Eb)));
	result = MR / sqrt(Ea * Eb);
	return result;
}

//****************************************************************************
extern "C" long double MRversusEavg(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double MR, result;
	long double Ea, Eb, pZa, pZb;
	if (vectorstoadd.size() == 2) {
		Ea = vec_part[vectorstoadd[0]].E();
		pZa = vec_part[vectorstoadd[0]].Pz();
		Eb = vec_part[vectorstoadd[1]].E();
		pZb = vec_part[vectorstoadd[1]].Pz();
	} else {
		cout << " betarazor can be evaluated only for two particles" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	MR =
			2
					* sqrt(
							(Ea * pZb - Eb * pZa) * (Ea * pZb - Eb * pZa)
									/ ((pZa - pZb) * (pZa - pZb)
											- (Ea - Eb) * (Ea - Eb)));
	result = 2 * MR / (Ea + Eb);
	return result;
}

//****************************************************************************
extern "C" long double Egeomavg(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double MR, result;
	long double Ea, Eb, pZa, pZb;
	if (vectorstoadd.size() == 2) {
		Ea = vec_part[vectorstoadd[0]].E();
		pZa = vec_part[vectorstoadd[0]].Pz();
		Eb = vec_part[vectorstoadd[1]].E();
		pZb = vec_part[vectorstoadd[1]].Pz();
	} else {
		cout << " Eratio can be evaluated only for two particles" << endl;
		exit(1);
	}
	//cout<<pt<<endl;

	result = sqrt(Ea * Eb);
	return result;
}

//****************************************************************************
extern "C" long double rapidity(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double rapidity;
	long double E, pZ;
	if (vectorstoadd.size() == 1) {
		E = vec_part[vectorstoadd[0]].E();
		pZ = vec_part[vectorstoadd[0]].Pz();
	} else {
		cout << " pT can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<pt<<endl;
	rapidity = 0.5 * log((E + pZ) / (E - pZ));
	return rapidity;
}

//****************************************************************************
extern "C" long double rapidityofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double rapidityofn, E, pZ;

	TLorentzVector sum;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	E = sum.E();
	pZ = sum.Pz();

	//cout<<pZ<<endl;
	rapidityofn = 0.5 * log((E + pZ) / (E - pZ));
	return sum.Rapidity();

}

//****************************************************************************
extern "C" long double xPDFmin(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, rootS, E, pZ, x1, x2;

	TLorentzVector sum;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	E = sum.E();
	pZ = sum.Pz();

	rootS = 7000.;

	x1 = (E + pZ) / rootS;
	x2 = (E - pZ) / rootS;

	res = min(x1, x2);

	return res;

}

//****************************************************************************
extern "C" long double xPDFmax(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, rootS, E, pZ, x1, x2;

	TLorentzVector sum;
	//puts("involved particles ");
	//printvector(vectorstoadd);

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	E = sum.E();
	pZ = sum.Pz();

	rootS = 7000.;

	x1 = (E + pZ) / rootS;
	x2 = (E - pZ) / rootS;

	res = max(x1, x2);

	return res;

}

//****************************************************************************
extern "C" long double E(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double E;
	if (vectorstoadd.size() == 1) {
		E = vec_part[vectorstoadd[0]].E();
		//cout << E << endl;

	} else {
		cout << " E can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout << E << endl;

	return E;
}

//****************************************************************************
extern "C" long double sizeofevent(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double number;
	long double ve[5];

	number = 0;
	cout << "fsdjkfjksdfhsdhfdsjk" << endl;

	cout << "size del vettore " << sizeof(vec_part) << endl;

	cout << "size del 0-esimo " << sizeof(vec_part[0]) << endl;
	cout << "size del 1-esimo " << sizeof(vec_part[1]) << endl;

	cout << "size del vettore " << sizeof(ve) << endl;

	cout << "size del 0-esimo " << sizeof(ve[0]) << endl;
	cout << "size del 1-esimo " << sizeof(ve[1]) << endl;

	//number=sizeof(vec_part)/sizeof(vec_part[0]);
	return number;
}

//****************************************************************************
extern "C" long double number(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	//cout << "sono dentro number" << endl;
	//printmatrixgen(matrixtoadd);
	long double res;
	if (matrixtoadd.empty() == false) {
		res = matrixtoadd.size();
		//cout << "the number is " << res << endl;
	} else if (matrixtoadd.empty() == true) {
		res = 0;
		//cout << "the number is " << res << endl;

	}
	//cout<< sizeof(&vec_part[0])<<endl;
	//number=sizeof(vec_part)/sizeof(vec_part[0]);
	//cout << res << endl;
	return res;
}

//****************************************************************************
extern "C" long double numberofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	//puts("sono dentro number");

	if (1 == matrixtoadd.size()) {

		long double number;
		vector<int> vectorstoadd = matrix2vector(matrixtoadd);

		number = vectorstoadd.size();
		return number;
	} else {
		cout
				<< "the number of components of a composite object can be computed onlty on one object at a time"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double numberratio(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	//puts("sono dentro number");

	if (2 == matrixtoadd.size()) {

		int number1;
		int number2;

		number1 = matrixtoadd[0].size();
		number2 = matrixtoadd[1].size();
		//		cout << number1 << endl;
		//		cout << number2 << endl;
		//		cout << exp(-fabs(log(number1/number2))) << endl;

		//		return exp(-fabs(log(number1 / number2)));
		long double ratio;
		ratio = (long double) min(number1, number2) / max(number1, number2);
		return ratio;
	} else {
		cout
				<< "the ratio of number of components of composites objects can be computed onlty on two objects at a time"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double centralpT15eta25(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 1) {
		long double eta, pT;
		eta = vec_part[vectorstoadd[0]].Eta();
		pT = vec_part[vectorstoadd[0]].Pt();

		if ((fabs(eta) < 2.5) && (pT > 15)) {
			return 1;

		} else {
			return 0;
		}

	} else {
		cout << "centrality can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double centralpT20eta25(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 1) {
		long double eta, pT;
		eta = vec_part[vectorstoadd[0]].Eta();
		pT = vec_part[vectorstoadd[0]].Pt();

		if ((fabs(eta) < 2.5) && (pT > 20)) {
			return 1;

		} else {
			return 0;
		}

	} else {
		cout << "centrality can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double centralpT20eta3(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 1) {
		long double eta, pT;
		eta = vec_part[vectorstoadd[0]].Eta();
		pT = vec_part[vectorstoadd[0]].Pt();

		if ((fabs(eta) < 3) && (pT > 20)) {
			return 1;

		} else {
			return 0;
		}

	} else {
		cout << "centrality can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double centralpT30eta3(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 1) {
		long double eta, pT;
		eta = vec_part[vectorstoadd[0]].Eta();
		pT = vec_part[vectorstoadd[0]].Pt();

		if ((fabs(eta) < 3) && (pT > 30)) {
			return 1;

		} else {
			return 0;
		}

	} else {
		cout << "centrality can be evaluated only on one particle" << endl;
		exit(1);
	}
}

//****************************************************************************
extern "C" long double eta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Eta();
	} else {
		cout << "eta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;
	return eta;
}

//****************************************************************************
extern "C" long double rapidity_FromROOT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Rapidity();
	} else {
		cout << "eta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;
	return eta;
}

//****************************************************************************
extern "C" long double yStar(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double yStar, eta1, eta2;

	if (vectorstoadd.size() == 2) {
		eta1 = vec_part[vectorstoadd[0]].Eta();
		eta2 = vec_part[vectorstoadd[1]].Eta();

	} else {
		cout << "chi can be evaluated only on a pair of particles" << endl;
		exit(1);
	}

	yStar = 0.5 * fabs(eta2 - eta1);

	return yStar;
}
//****************************************************************************
extern "C" long double yB(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double yB, eta1, eta2;

	if (vectorstoadd.size() == 2) {
		eta1 = vec_part[vectorstoadd[0]].Eta();
		eta2 = vec_part[vectorstoadd[1]].Eta();

	} else {
		cout << "chi can be evaluated only on a pair of particles" << endl;
		exit(1);
	}

	yB = 0.5 * fabs(eta2 + eta1);

	return yB;
}

//****************************************************************************
extern "C" long double chi_eta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double chi, eta1, eta2;

	if (vectorstoadd.size() == 2) {
		eta1 = vec_part[vectorstoadd[0]].Eta();
		eta2 = vec_part[vectorstoadd[1]].Eta();

	} else {
		cout << "chi can be evaluated only on a pair of particles" << endl;
		exit(1);
	}

	chi = exp(fabs(eta2 - eta1));

	return chi;
}

//****************************************************************************
extern "C" long double tagjet(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double yesorno;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta, pT, energy;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

	} else {
		cout << "eta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;

	if (energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5) {
		yesorno = 1;
	} else {
		yesorno = 0;
	}

	return yesorno;
}

//****************************************************************************
extern "C" long double eta_nottagjet(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double result;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta, pT, energy;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

	} else {
		cout << "taggines can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;

	if (energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5) {
		result = -100;
	} else {
		result = eta;
	}

	return result;
}

//****************************************************************************
extern "C" long double PTnottagjet(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double pt;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta, pT, energy;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

	} else {
		cout << "eta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;

	if (energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5) {
		pt = -100;
	} else {
		pt = pT;
	}

	return pt;
}

//****************************************************************************
extern "C" long double PTtagjet(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);
	long double pt;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta, pT, energy;

	if (vectorstoadd.size() == 1) {
		eta = vec_part[vectorstoadd[0]].Eta();
		energy = vec_part[vectorstoadd[0]].E();
		pT = vec_part[vectorstoadd[0]].Pt();

	} else {
		cout << "eta can be evaluated only on one particle" << endl;
		exit(1);
	}
	//cout<<eta<<endl;

	if (energy > 800 && pT > 40 && fabs(eta) > 3 && fabs(eta) < 5) {
		pt = pT;
	} else {
		pt = -100;
	}

	return pt;
}

//****************************************************************************
extern "C" long double pTofLightParton(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double pT;
	int vect = 0;
	TLorentzVector sum, visible;
	vector<int> lightpartons;
	lightpartons.push_back(21);
	lightpartons.push_back(1);
	lightpartons.push_back(2);
	lightpartons.push_back(3);
	lightpartons.push_back(-4);
	lightpartons.push_back(-1);
	lightpartons.push_back(-2);
	lightpartons.push_back(-3);
	lightpartons.push_back(-4);

	//cout << vect << " " << part[vect]->PID << endl;

	//	float length;
	//	length=sizeof(vec_part[8]);
	//	cout<<length<<endl;
	//  int pid;
	//  pid = part[0]->PID;
	//  cout << pid;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		//warning MG showered files contains some status 2 particle, which is not present in pythia output or in fastjet output

		if ((part[vectorstoadd[vect]]->Status > 0)
				&& (part[vectorstoadd[vect]]->Status != 2)) {

			if (isinclass(lightpartons, part[vectorstoadd[vect]]->PID)
					== true) {

				// it is a parton, thus compute the pT
				pT = vec_part[vectorstoadd[vect]].Pt();
			} else if (isinclass(lightpartons, part[vectorstoadd[vect]]->PID)
					== false) {

				// it is NOT a parton, thus give unphysical -100
				//pT = -100;

				// it is NOT a parton, thus do not do anything

			}
		}

	}

	//cout<<pT<<endl;
	return pT;
}

//****************************************************************************
extern "C" long double etaofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double etaofn;
	int vect = 0;
	TLorentzVector sum;

	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	etaofn = sum.Eta();
	//cout<<pTofn<<endl;
	return etaofn;
}
//****************************************************************************
extern "C" long double CDphi(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	long double deltaphi;

	if (vectorstoadd.size() == 2) {
		deltaphi = vec_part[vectorstoadd[0]].Vect().DeltaPhi(
				vec_part[vectorstoadd[1]].Vect()); // get deltaPhi angle between v1 and v2
	} else {
		cout << "deltaphi can be evaluated only for two fourvectors" << endl;
		exit(1);
	}

	//cout<<eta<<endl;
	res = cos(deltaphi);
	return res;
}
//****************************************************************************
extern "C" long double CDphiChDeta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta1, eta2, deta, res;

	if (vectorstoadd.size() == 2) {
		eta1 = vec_part[vectorstoadd[0]].Eta();
		eta2 = vec_part[vectorstoadd[1]].Eta();
		deta = fabs(eta1 - eta2);
	} else {
		cout
				<< "cos(deltaphi) - cosh(delta_eta) can be evaluated only on two particles"
				<< endl;
		exit(1);
	}

	long double deltaphi;

	if (vectorstoadd.size() == 2) {
		deltaphi = vec_part[vectorstoadd[0]].Vect().DeltaPhi(
				vec_part[vectorstoadd[1]].Vect()); // get angle between v1 and v2
	} else {
		cout << "deltaphi can be evaluated only for two fourvectors" << endl;
		exit(1);
	}

	//cout<<eta<<endl;
	res = -cos(deltaphi) + cosh(deta);
	return res;
}

//****************************************************************************
extern "C" long double delta_rho(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	long double res;

	res = CDphiChDeta(matrixtoadd, vec_part, part);

	return res;
}

//****************************************************************************
extern "C" long double delta_eta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double eta1, eta2, res;

	if (vectorstoadd.size() == 2) {
		eta1 = vec_part[vectorstoadd[0]].Eta();
		eta2 = vec_part[vectorstoadd[1]].Eta();
		res = fabs(eta1 - eta2);
	} else {
		cout << "delta_eta can be evaluated only on two particles" << endl;
		exit(1);
	}
	//cout<<eta<<endl;
	return res;
}

//****************************************************************************
extern "C" long double delta_y_t_minus_tbar(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double etatop, etatopbar, res;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;

		if (pid1 == 6 || pid1 == 66) {
			etatop = vec_part[vectorstoadd[0]].Rapidity();
		} else if (pid1 == -6 || pid1 == -66) {
			etatopbar = vec_part[vectorstoadd[0]].Rapidity();
		}
		if (pid2 == 6 || pid2 == 66) {
			etatop = vec_part[vectorstoadd[1]].Rapidity();
		} else if (pid2 == -6 || pid2 == -66) {
			etatopbar = vec_part[vectorstoadd[1]].Rapidity();
		}

		res = etatop - etatopbar;

	} else {
		cout << "delta_y can be evaluated only on two particles" << endl;
		exit(1);
	}
	//cout<<eta<<endl;
	return res;
}

//****************************************************************************
extern "C" long double delta_eta_t_minus_tbar(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double etatop, etatopbar, res;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;

		if (pid1 == 6) {
			etatop = vec_part[vectorstoadd[0]].Eta();
		} else if (pid1 == -6) {
			etatopbar = vec_part[vectorstoadd[0]].Eta();
		}
		if (pid2 == 6) {
			etatop = vec_part[vectorstoadd[1]].Eta();
		} else if (pid2 == -6) {
			etatopbar = vec_part[vectorstoadd[1]].Eta();
		}

		res = etatop - etatopbar;

	} else {
		cout << "delta_y can be evaluated only on two particles" << endl;
		exit(1);
	}
	//cout<<eta<<endl;
	return res;
}

//****************************************************************************
extern "C" long double pseudorapidity(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double costheta;

	if (vectorstoadd.size() == 1) {
		costheta = vec_part[vectorstoadd[0]].Vect().CosTheta();
	} else {
		cout << "costheta can be evaluated only on one particle" << endl;
		exit(1);
	}

	long double eta = 0.5 * log((1 + costheta) / (1 - costheta)); //} else {
	return eta;
	//}
}

//****************************************************************************
extern "C" long double costheta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double costheta;

	if (vectorstoadd.size() == 1) {
		costheta = vec_part[vectorstoadd[0]].Vect().CosTheta();
	} else {
		cout << "costheta can be evaluated only on one particle" << endl;
		exit(1);
	}

	return costheta; //} else {

	//}
}

//****************************************************************************
extern "C" long double sintheta(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double sintheta;

	if (vectorstoadd.size() == 1) {
		sintheta = sin(vec_part[vectorstoadd[0]].Vect().Theta());
	} else {
		cout << "costheta can be evaluated only on one particle" << endl;
		exit(1);
	}

	return sintheta; //} else {

	//}
}

//****************************************************************************
extern "C" long double cos_delta_phi(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double cosdeltaphi, pxA, pxB, pyA, pyB;

	if (vectorstoadd.size() == 2) {

		pxA = vec_part[vectorstoadd[0]].Px();
		pyA = vec_part[vectorstoadd[0]].Py();

		pxB = vec_part[vectorstoadd[1]].Px();
		pyB = vec_part[vectorstoadd[1]].Py();

		cosdeltaphi = (pxA * pxB + pyA * pyB) / sqrt(pxA * pxA + pyA * pyA)
				/ sqrt(pxB * pxB + pyB * pyB);
		//fabs(
		//	vec_part[vectorstoadd[0]].Vect().DeltaPhi(
		//		vec_part[vectorstoadd[1]].Vect())); // get angle between v1 and v2
	} else {
		cout << "deltaphi can be evaluated only for two fourvectors" << endl;
		exit(1);
	}

	return cosdeltaphi; //} else {

	//}
}

//****************************************************************************
extern "C" long double delta_phi(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double deltaphi;

	if (vectorstoadd.size() == 2) {
		deltaphi = fabs(
				vec_part[vectorstoadd[0]].Vect().DeltaPhi(
						vec_part[vectorstoadd[1]].Vect())); // get angle between v1 and v2
	} else {
		cout << "deltaphi can be evaluated only for two fourvectors" << endl;
		exit(1);
	}

	return deltaphi; //} else {

	//}
}

//****************************************************************************
extern "C" long double deltaR(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double deltaR;

	if (vectorstoadd.size() == 2) {
		deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(
				vec_part[vectorstoadd[1]].Vect());
	} else {
		cout << "deltaR can be evaluate only for 2 fourvectors" << endl;
		exit(1);
	}

	return deltaR; //} else {

	//}
}

//****************************************************************************
extern "C" long double deltaRjetnonjetUNDER(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double deltaR;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		//cout << pid1 << " " << pid2 << endl;

		if ((is_a_parton(pid1) == false) && (is_a_parton(pid2) == true)
				|| (is_a_parton(pid2) == false) && (is_a_parton(pid1) == true)) {

			deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(
					vec_part[vectorstoadd[1]].Vect());

		} else {
			// throw and output out of the physical range as garbage container.
			// as this observable is used mostly for isolation cuts it make sense to have a garbage bin
			// in a range typically accepted by the isolation cuts
			// the max \Delta\R  is -100 which is a garbage bin is
			deltaR = -100;
		}
	} else {
		cout << "deltaR can be evaluate only for 2 fourvectors" << endl;
		exit(1);
	}

	return deltaR; //} else {

	//}
}

//****************************************************************************
extern "C" long double deltaRjetX(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double deltaR;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		//cout << pid1 << " " << pid2 << endl;

		if (((pid2 == 93) && (pid1 != 93)) || ((pid1 == 93) && (pid2 != 93))) {

			deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(
					vec_part[vectorstoadd[1]].Vect());

		} else {
			// throw and output out of the physical range as garbage container.
			// as this observable is used mostly for isolation cuts it make sense to have a garbage bin
			// in a range typically accepted by the isolation cuts
			// the max \Delta\eta^2 is 25 so a good garbage bin is
			deltaR = 35;
		}
	} else {
		cout << "deltaR can be evaluate only for 2 fourvectors" << endl;
		exit(1);
	}

	return deltaR; //} else {

	//}
}

//****************************************************************************
extern "C" long double deltaRdiff(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the total invariant mass of N particles
	// particle involved are those in the event with IDs listed in  the vector vectorstoadd
	//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double deltaR;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		//cout << pid1 << " " << pid2 << endl;

		if (pid1 != pid2) {

			deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(
					vec_part[vectorstoadd[1]].Vect());
		} else {
			//throw and output out of the physical range as garbage container
			deltaR = -1;
		}
	} else {
		cout << "deltaR can be evaluate only for 2 fourvectors" << endl;
		exit(1);
	}

	return deltaR; //} else {

	//}
}

//****************************************************************************
extern "C" long double mT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the angle of the cosine brtween two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double cTij, deltaphi, mT, Ei, Ej;

	if (vectorstoadd.size() == 2) {

		deltaphi = vec_part[vectorstoadd[0]].Vect().DeltaPhi(
				vec_part[vectorstoadd[1]].Vect());

	} else {
		cout << "mT can be evaluate only for two fourvectors" << endl;
		exit(1);
	}
	cTij = cos(deltaphi);
	Ei = vec_part[vectorstoadd[0]].E();
	Ej = vec_part[vectorstoadd[1]].E();
	mT = sqrt(2 * Ei * Ej * (1 - cTij));
	return mT; //} else {

	//}
}

extern "C" long double mTmassless(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the angle of the cosine brtween two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double cTij, deltaphi, mT, Pti, Ptj;

	if (vectorstoadd.size() == 2) {

		deltaphi = vec_part[vectorstoadd[0]].Vect().DeltaPhi(
				vec_part[vectorstoadd[1]].Vect());

	} else {
		cout << "mT can be evaluate only for two fourvectors" << endl;
		exit(1);
	}
	cTij = cos(deltaphi);
	Pti = vec_part[vectorstoadd[0]].Pt();
	Ptj = vec_part[vectorstoadd[1]].Pt();
	mT = sqrt(2 * Pti * Ptj * (1 - cTij));
	return mT; //} else {

	//}
}

//****************************************************************************
extern "C" long double mT_blind(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************
	//***** PIDs must contain the neutrinos AND the particle to be considered as visible decay products
	//***** then sum al the visibile particles into one vector and
	//***** compute the transverse mass assuming that there is a unique mET source, though there could be many neutrinos

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxV, pyV;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	// for Delphes LHCO converted to LHE
	neutrinos.push_back(121416);

	// find the Px and Py of the system of invisible particles

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// find the Px and Py of the system of invisible particles
	// visible particles are those not being added to the list of neutrinos above

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == false) {
			pxV = vec_part[vectorstoadd[kpart]].Px();
			pyV = vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// ****** Compute mT *******
	// as in the ATLAS TDR pg 706 because I want to reproduce their plot
	float abs_pTV, abs_pTmiss, pTVdotpTmiss;

	abs_pTV = sqrt(pxV * pxV + pyV * pyV);
	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);
	pTVdotpTmiss = pxmiss * pxV + pymiss * pyV;
	res = sqrt(2 * (abs_pTmiss * abs_pTV - pTVdotpTmiss));

	return res;

}

//****************************************************************************
extern "C" long double mT_HDO(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************
	//***** PIDs must contain the neutrinos AND the particle to be considered as visible decay products
	//***** then sum al the visibile particles into one vector and
	//***** compute the transverse mass assuming that there is a unique mET source, though there could be many neutrinos

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxV, pyV;
	pxmiss = 0;
	pymiss = 0;
	pxV = 0;
	pyV = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	neutrinos.push_back(121416);

	// find the Px and Py of the system of invisible particles

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// find the Px and Py of the system of invisible particles
	// visible particles are those not being added to the list of neutrinos above

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == false) {
			cout << pids[kpart] << endl;
			pxV = pxV + vec_part[vectorstoadd[kpart]].Px();
			//cout<< pids[kpart] << << pxV endl;
			pyV = pyV + vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// ****** Compute mT *******
	// as in the ATLAS TDR pg 706 because I want to reproduce their plot
	float abs_pTV, abs_pTmiss, pTVdotpTmiss;

	abs_pTV = sqrt(pxV * pxV + pyV * pyV);
	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);
	pTVdotpTmiss = pxmiss * pxV + pymiss * pyV;
	res = sqrt(2 * (abs_pTmiss * abs_pTV - pTVdotpTmiss));

	cout << "pxV:" << pxV << "pyV:" << pyV << " pxmiss:" << pxmiss << " pymiss:"
			<< pymiss << " res:" << res << endl;

	return res;

}

//****************************************************************************
extern "C" long double s_min_basic(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************
	//***** PIDs must contain the neutrinos AND the particle to be considered as visible decay products
	//***** then sum al the visibile particles into one vector and
	//***** compute the transverse mass assuming that there is a unique mET source, though there could be many neutrinos

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxV, pyV, pzV, eV;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	neutrinos.push_back(121416);

	// find the Px and Py of the system of invisible particles

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// find the Px and Py of the system of invisible particles
	// visible particles are those not being added to the list of neutrinos above

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == false) {
			eV = vec_part[vectorstoadd[kpart]].E();
			pxV = vec_part[vectorstoadd[kpart]].Px();
			pyV = vec_part[vectorstoadd[kpart]].Py();
			pzV = vec_part[vectorstoadd[kpart]].Pz();
		}
	}

	// ****** Compute \sqrt{s}_min as 1006.0653 eq.(1.2) for MissingMass=0 *******

	float abs_pTV, abs_pTmiss, pTVdotpTmiss;

	abs_pTV = sqrt(pxV * pxV + pyV * pyV);
	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);

	res = sqrt(eV * eV - pzV * pzV) + abs_pTmiss;
	return res;

}

//****************************************************************************
extern "C" long double mT_approx(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************
	//***** PIDs must contain the neutrinos AND the particle to be considered as visible decay products
	//***** then sum al the visibile particles into one vector and
	//***** compute the transverse mass assuming that there is a unique mET source, though there could be many neutrinos

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;
	TLorentzVector pll;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxV, pyV;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	// for Delphes LHCO converted to LHE
	neutrinos.push_back(121416);

	// find the Px and Py of the system of invisible particles

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// find the Px and Py of the system of invisible particles
	// visible particles are those not being added to the list of neutrinos above

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == false) {

			pll = pll + vec_part[vectorstoadd[kpart]];

		}
	}

	// ****** Compute mT *******
	// as in hep-ph/0402254 because I want to reproduce their plot
	float abs_pTV, abs_pTmiss, pTVdotpTmiss, mll, ETll, ETnunu, absQ_pTllnunu;
	mll = sqrt(pll * pll);
	abs_pTV = fabs(pll.Pt());
	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);

	ETll = sqrt(abs_pTV * abs_pTV + mll * mll);
	ETnunu = sqrt(abs_pTmiss * abs_pTmiss + mll * mll);

	absQ_pTllnunu = (pll.Px() + pxmiss) * (pll.Px() + pxmiss)
			+ (pll.Py() + pymiss) * (pll.Py() + pymiss);

	res = sqrt((ETll + ETnunu) * (ETll + ETnunu) - absQ_pTllnunu);

	return res;

}

//****************************************************************************
extern "C" long double mT_approx_smeared(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************
	//***** PIDs must contain the neutrinos AND the particle to be considered as visible decay products
	//***** then sum al the visibile particles into one vector and
	//***** compute the transverse mass assuming that there is a unique mET source, though there could be many neutrinos

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, smear, sigma, HT;
	HT = 0;
	TLorentzVector pll;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxV, pyV;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	// for Delphes LHCO converted to LHE
	neutrinos.push_back(121416);

	// find the Px and Py of the system of invisible particles

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// find the Px and Py of the system of invisible particles
	// visible particles are those not being added to the list of neutrinos above

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == false) {

			pll = pll + vec_part[vectorstoadd[kpart]];
			HT = HT + fabs(vec_part[vectorstoadd[kpart]].Pt());
		}
	}

	// ****** Compute mT *******
	// as in hep-ph/0402254 because I want to reproduce their plot
	float abs_pTV, abs_pTmiss, pTVdotpTmiss, mll, ETll, ETnunu, absQ_pTllnunu,
			pxmissbuffer, pymissbuffer;
	mll = sqrt(pll * pll);
	abs_pTV = fabs(pll.Pt());

	// *******************************
	sigma = 1 * 0.5 * sqrt(HT);
	// *******************************
	TRandom *r3 = new TRandom3(0);

	smear = r3->Gaus(0, sigma);
	//cout<<smear<<"   ";
	pxmissbuffer = pxmiss + smear;
	smear = r3->Gaus(0, sigma);
	//cout<<smear<<endl;
	pymissbuffer = pymiss + smear;

	pxmiss = pxmissbuffer;
	pymiss = pymissbuffer;

	abs_pTmiss = sqrt(pxmiss * pxmiss + pymiss * pymiss);

	ETll = sqrt(abs_pTV * abs_pTV + mll * mll);
	ETnunu = sqrt(abs_pTmiss * abs_pTmiss + mll * mll);

	absQ_pTllnunu = (pll.Px() + pxmiss) * (pll.Px() + pxmiss)
			+ (pll.Py() + pymiss) * (pll.Py() + pymiss);

	res = sqrt((ETll + ETnunu) * (ETll + ETnunu) - absQ_pTllnunu);

	return res;

}

//
//****************************************************************************
extern "C" long double mT2_2nu_ideal(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res;

	//pa, pb = {mass, px, py}
	//pmiss  = {NULL, pxmiss, pymiss}
	//mn     = invisible particle mass

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pxA, pxB, pyA, pyB;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;
	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	// for Delphes LHCO converted to LHE
	neutrinos.push_back(121416);

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// take the two leptons and assing A and B their x and y momementa
	// pa, pb = {mass, px, py}
	// by definition A is the positve PID lepton and B has negative PID

	vector<int> leptonsA;
	vector<int> leptonsB;
	leptonsA.push_back(11);
	leptonsA.push_back(13);
	leptonsA.push_back(15);
	leptonsB.push_back(-11);
	leptonsB.push_back(-13);
	leptonsB.push_back(-15);

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(leptonsA, pids[kpart]) == true) {
			pxA = vec_part[vectorstoadd[kpart]].Px();
			pyA = vec_part[vectorstoadd[kpart]].Py();
		}
		if (isinclass(leptonsB, pids[kpart]) == true) {
			pxB = vec_part[vectorstoadd[kpart]].Px();
			pyB = vec_part[vectorstoadd[kpart]].Py();
		}
	}

	// ****** Compute mT2 *******

	double pa[3] = { 0, pxA, pyA };
	double pb[3] = { 0, pxB, pyB };
	double pmiss[3] = { 0, pxmiss, pymiss };
	double mn = 0.;

	mt2_bisect::mt2 mt2_event;

	mt2_event.set_momenta(pa, pb, pmiss);
	mt2_event.set_mn(mn);
	//	mt2_event.print();

	res = mt2_event.get_mt2();

	return res;

}

//****************************************************************************
extern "C" long double eta_aa_jjAVG(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, eta1, eta2, etaGammaGamma;
	vector<TLorentzVector> jets, gammas;
	TLorentzVector j1, j2, gamma1, gamma2, sum;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	vector<int> jetPID, gammasPID;
	vector<float> pTsJ, pTsGamma;

	jetPID.push_back(55);
	jetPID.push_back(-55);
	jetPID.push_back(211234);

	gammasPID.push_back(22);

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(jetPID, pids[kpart]) == true) {

			pTsJ.push_back(vec_part[vectorstoadd[kpart]].Pt());
			jets.push_back(vec_part[vectorstoadd[kpart]]);

		}
	}

	if (pTsJ.size() > 1) {

		vector<int> idjets;
		idjets = MaxKElementsInAVector(pTsJ, 2);

		//PrintVector(std::cout, idjets, " ");
		//cout << endl;
		//	PrintVector(std::cout, pTsJ, " ");
		//	cout << endl;

		j1 = jets[idjets[0]];
		j2 = jets[idjets[1]];

//		cout << "The jet pT1 and pT2 are " << j1.Pt() << " " << j2.Pt() << endl;

		eta1 = j1.Eta();
		eta2 = j2.Eta();

	} else {
		eta1 = 100;
		eta2 = 0;
	}

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(gammasPID, pids[kpart]) == true) {

			pTsGamma.push_back(vec_part[vectorstoadd[kpart]].Pt());
			gammas.push_back(vec_part[vectorstoadd[kpart]]);

		}
	}

	if (pTsGamma.size() > 1) {
		vector<int> idgammas;
		idgammas = MaxKElementsInAVector(pTsGamma, 2);

		//	PrintVector(std::cout, idgammas, " ");
		//		cout << endl;
		//	PrintVector(std::cout, pTsGamma, " ");
		//cout << endl;

		gamma1 = gammas[idgammas[0]];
		gamma2 = gammas[idgammas[1]];

//		cout << "The Gamma pT1 and pT2 are " << gamma1.Pt() << " " << gamma2.Pt() << endl;

		sum = gamma2 + gamma1;

		etaGammaGamma = sum.Eta();
	} else {
		etaGammaGamma = 200;
	}
	res = fabs(0.5 * (eta1 + eta2) - etaGammaGamma);

	return res;

}

//****************************************************************************
extern "C" long double DeltaPhi_aa_jj(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************

	bool debug;
	debug = true;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, eta1, eta2, etaGammaGamma;
	vector<TLorentzVector> jets, gammas;
	TLorentzVector j1, j2, gamma1, gamma2, sumjj, sumgammagamma;

	vector<TRootLHEFParticle*> jetLHE, gammaLHE;

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

		if (debug) {

			cout << part[vectorstoadd[kpart]]->PID << " "
					<< vec_part[vectorstoadd[kpart]].Pt() << " "
					<< vec_part[vectorstoadd[kpart]].Px() << " "
					<< vec_part[vectorstoadd[kpart]].Py() << endl;
		}

	}

	vector<int> jetPID, gammasPID;
	vector<float> pTsJ, pTsGamma;

	jetPID.push_back(55);
	jetPID.push_back(-55);
	jetPID.push_back(211234);

	gammasPID.push_back(22);

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(jetPID, pids[kpart]) == true) {

			pTsJ.push_back(vec_part[vectorstoadd[kpart]].Pt());
			jets.push_back(vec_part[vectorstoadd[kpart]]);
			jetLHE.push_back(part[vectorstoadd[kpart]]);

		}
	}

	if (pTsJ.size() > 1) {

		vector<int> idjets;
		idjets = MaxKElementsInAVector(pTsJ, 2);

		if (debug) {
			PrintVector(std::cout, idjets, " ");
			cout << endl;
			PrintVector(std::cout, pTsJ, " ");
			cout << endl;
		}

		j1 = jets[idjets[0]];
		j2 = jets[idjets[1]];

		sumjj = j1 + j2;

		if (debug) {
			cout << "The jet pT1 and pT2 are " << j1.Pt() << " " << j2.Pt()
					<< endl;
			cout << "The PIDs of these particles are " << jetLHE[idjets[0]]->PID
					<< " " << jetLHE[idjets[1]]->PID << endl;
		}
	} else {
		res = -20;
	}

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(gammasPID, pids[kpart]) == true) {

			pTsGamma.push_back(vec_part[vectorstoadd[kpart]].Pt());
			gammas.push_back(vec_part[vectorstoadd[kpart]]);
			gammaLHE.push_back(part[vectorstoadd[kpart]]);

		}
	}

	if (pTsGamma.size() > 1) {
		vector<int> idgammas;
		idgammas = MaxKElementsInAVector(pTsGamma, 2);

		if (debug) {
			PrintVector(std::cout, idgammas, " ");
			cout << endl;
			PrintVector(std::cout, pTsGamma, " ");
			cout << endl;
		}

		gamma1 = gammas[idgammas[0]];
		gamma2 = gammas[idgammas[1]];

		if (debug) {
			cout << "The Gamma pT1 and pT2 are " << gamma1.Pt() << " "
					<< gamma2.Pt() << endl;

			cout << "The PIDs of these particles are "
					<< gammaLHE[idgammas[0]]->PID << " "
					<< gammaLHE[idgammas[1]]->PID << endl;

		}

		sumgammagamma = gamma2 + gamma1;

	} else {
		res = -20;
	}

	res = fabs(sumgammagamma.Vect().DeltaPhi(sumjj.Vect()));

	if (debug) {
		float dphij1j2, dphigamma1gamma2, dphij1gamma2, dphij1gamma1,
				dphij2gamma1, dphij2gamma2;

		dphij1j2 = fabs(j1.Vect().DeltaPhi(j2.Vect())); // get angle between v1 and v2

		dphigamma1gamma2 = fabs(gamma1.Vect().DeltaPhi(gamma2.Vect()));

		dphij1gamma1 = fabs(j1.Vect().DeltaPhi(gamma1.Vect()));

		dphij1gamma2 = fabs(j1.Vect().DeltaPhi(gamma2.Vect()));

		dphij2gamma1 = fabs(j2.Vect().DeltaPhi(gamma1.Vect()));

		dphij2gamma2 = fabs(j2.Vect().DeltaPhi(gamma2.Vect()));

		cout << "j1j2: " << dphij1j2 << ", j1a2: " << dphij1gamma2 << ", j1a1: "
				<< dphij1gamma1 << ", j2a1: " << dphij2gamma1 << ", j2a2: "
				<< dphij2gamma2 << ", a1a2: " << dphigamma1gamma2 << endl;
	}
//	res = fabs(sumjj.Eta() - sumgammagamma.Eta());

	return res;

}

//****************************************************************************
extern "C" long double mT_ll_mET(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {

	//****************************************************************************

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double res, pTll, mET, deltaphi;
	TLorentzVector mETfv, llfv;

	//pa, pb = {mass, px, py}
	//pmiss  = {NULL, pxmiss, pymiss}
	//mn     = invisible particle mass

	vector<int> pids;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		pids.push_back(part[vectorstoadd[kpart]]->PID);

	}

	double pxmiss, pymiss, pXll, pYll;
	pxmiss = 0;
	pymiss = 0;

	vector<int> neutrinos;

	neutrinos.push_back(12);
	neutrinos.push_back(14);
	neutrinos.push_back(16);
	neutrinos.push_back(-12);
	neutrinos.push_back(-14);
	neutrinos.push_back(-16);
	// for Delphes LHCO converted to LHE
	neutrinos.push_back(121416);

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(neutrinos, pids[kpart]) == true) {
			pxmiss += vec_part[vectorstoadd[kpart]].Px();
			pymiss += vec_part[vectorstoadd[kpart]].Py();

			mETfv = mETfv + vec_part[vectorstoadd[kpart]];

		}
	}

	// take the two leptons and assing A and B their x and y momementa
	// pa, pb = {mass, px, py}
	// by definition A is the positve PID lepton and B has negative PID

	vector<int> leptons;
	leptons.push_back(11);
	leptons.push_back(13);
	leptons.push_back(15);
	leptons.push_back(-11);
	leptons.push_back(-13);
	leptons.push_back(-15);

	pXll = 0;
	pYll = 0;
	int lepton_counter = 0;

	for (int kpart = 0; kpart < vectorstoadd.size(); kpart++) {

		if (isinclass(leptons, pids[kpart]) == true) {
			pXll = vec_part[vectorstoadd[kpart]].Px();
			pYll = vec_part[vectorstoadd[kpart]].Py();
			lepton_counter++;

			llfv = llfv + vec_part[vectorstoadd[kpart]];

		}

	}

	// ****** Compute mT *******

	pTll = llfv.Pt();
	mET = mETfv.Pt();

	deltaphi = fabs(llfv.Vect().DeltaPhi(mETfv.Vect())); // get angle between v1 and v2

	res = sqrt(2. * mET * pTll * (1. - cos(deltaphi)));

	return res;

}

//****************************************************************************
extern "C" long double lifetime(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the cosine of the  angle between two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double tau;

	if (vectorstoadd.size() == 1) {
		tau = part[vectorstoadd[0]]->LifeTime;
	} else {
		cout << "lifetime can be read only for 1 single paticle" << endl;
		exit(1);
	}

	return tau; //} else {

	//}
}

//****************************************************************************
extern "C" long double helicity(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the cosine of the  angle between two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double hel;

	if (vectorstoadd.size() == 1) {
		hel = part[vectorstoadd[0]]->Spin;
	} else if (vectorstoadd.size() == 2) {
		hel = pow(2, part[vectorstoadd[0]]->Spin)
				+ pow(2, part[vectorstoadd[1]]->Spin);
	} else {
		cout << "helicity can be read only for 1 single paticle" << endl;
		exit(1);
	}

	return hel; //} else {

	//}
}

//****************************************************************************
extern "C" long double cosij(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
	//****************************************************************************
	// compute the cosine of the  angle between two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double cij, a;

	if (vectorstoadd.size() == 2) {
		a = vec_part[vectorstoadd[0]].Vect().Angle(
				vec_part[vectorstoadd[1]].Vect()); // get angle between v1 and v2
	} else {
		cout << "cosij can be evaluate only for two fourvectors" << endl;
		exit(1);
	}
	cij = cos(a);
	return cij; //} else {

	//}
}

//****************************************************************************
extern "C" long double cosijMotherChildWPlus(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector) {
	//****************************************************************************
	// compute the cosine of the  angle between the sum of two particles (the mother) and the child with PID>0, the fermion
	// this function assumes that only one pair of children exists. i.e. WW has to decay lnujj and not jj jj beacause no check on the charge of the pair is made.

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	//	cout << "******************************" << vectorstoadd.size()<<endl;
	//	printvector(vectorstoadd);

	TLorentzVector sum_of_childrenChargePlus, PIDplus;

	int nmax = vectorstoadd.size();

	//---->loop over *particles*<------
	for (Int_t jcount = 0; jcount < nmax; ++jcount) {

		int status;
		status = partvector[vectorstoadd[jcount]]->Status;
		if (status == -23) {

			float charge = pid23timescharge(
					partvector[vectorstoadd[jcount]]->PID) / 3.0;
			//			cout<<partvector[vectorstoadd[jcount]]->PID<<" "<<charge<<endl;

			if (charge > 0) {
				//				cout << "una particella positiva e` " << vectorstoadd[jcount] << endl;
				sum_of_childrenChargePlus = sum_of_childrenChargePlus
						+ vec_partvector[vectorstoadd[jcount]];

				if (partvector[vectorstoadd[jcount]]->PID > 0) {
					//					cout << "il fermione carico positivo e` " << vectorstoadd[jcount] << endl;
					PIDplus = vec_partvector[vectorstoadd[jcount]];
				}

			}

		}
	}

	TVector3 bCoM;
	TLorentzVector motherLAB;
	motherLAB = sum_of_childrenChargePlus;

	bCoM = -sum_of_childrenChargePlus.BoostVector(); // b=(x/t,y/t,z/t)

	PIDplus.Boost(bCoM);

	long double cij, a;

	a = PIDplus.Vect().Angle(motherLAB.Vect()); // get angle between v1 and v2

	cij = cos(a);

	//	cout << "angolo e coseno fermione-W "<< a << " " << cij << endl;
	return cij;

}

//****************************************************************************
extern "C" long double cosMotherChildbQuarkInMotherRestFrame(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector) {
	//****************************************************************************
	// compute the cosine of the  angle between the sum of two particles (the mother) and the child with PID>0, the fermion
	// this function assumes that only one pair of children exists. i.e. WW has to decay lnujj and not jj jj beacause no check on the charge of the pair is made.

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	//	cout << "******************************" << vectorstoadd.size()<<endl;
	//	printvector(vectorstoadd);

	TLorentzVector sum_of_children, bquarkmomentum;

	int nmax = vectorstoadd.size();

	if (vectorstoadd.size() == 2) {
		//---->loop over *particles*<------
		for (Int_t jcount = 0; jcount < nmax; ++jcount) {

			sum_of_children = sum_of_children
					+ vec_partvector[vectorstoadd[jcount]];

			if (abs(partvector[vectorstoadd[jcount]]->PID) != 24) {
				int outputPID = abs(partvector[vectorstoadd[jcount]]->PID);
				//cout << "il fermione (il non W) e` " << outputPID << endl;
				bquarkmomentum = vec_partvector[vectorstoadd[jcount]];
			}

		}

		TVector3 bCoM;
		TLorentzVector motherLAB;
		motherLAB = sum_of_children;

		bCoM = -sum_of_children.BoostVector(); // b=(x/t,y/t,z/t)

		bquarkmomentum.Boost(bCoM);

		long double cij, a;

		a = bquarkmomentum.Vect().Angle(motherLAB.Vect()); // get angle between v1 and v2

		cij = cos(a);

		//	cout << "angolo e coseno fermione-W "<< a << " " << cij << endl;
		return cij;
	} else if (vectorstoadd.size() != 2) {

		cout << "the cosine in the CoM can be computed only for two particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double cosMotherChildbQuarkInMotherRestFrameOfN(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_partvector,
		vector<TRootLHEFParticle*> partvector) {
	//****************************************************************************
	// compute the cosine of the  angle between the sum of N particles (the mother) and the child with PID=5 or PID=-5
	//****************************************************************************

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	//	cout << "******************************" << vectorstoadd.size()<<endl;
	//	printvector(vectorstoadd);

	TLorentzVector sum_of_children, bquarkmomentum;

	int nmax = vectorstoadd.size();

	if (vectorstoadd.size() > 1) {
		//---->loop over *particles*<------
		for (Int_t jcount = 0; jcount < nmax; ++jcount) {

			sum_of_children = sum_of_children
					+ vec_partvector[vectorstoadd[jcount]];

			if (abs(partvector[vectorstoadd[jcount]]->PID) == 5
					|| (partvector[vectorstoadd[jcount]]->PID == 13)
					|| (partvector[vectorstoadd[jcount]]->PID == -1000006)) {
				int outputPID = partvector[vectorstoadd[jcount]]->PID;
				//cout << "il fermione (il non-W non-g) e` " << outputPID << endl;
				bquarkmomentum = vec_partvector[vectorstoadd[jcount]];
			}

		}

		TVector3 bCoM;
		TLorentzVector motherLAB;
		motherLAB = sum_of_children;

		bCoM = -sum_of_children.BoostVector(); // b=(x/t,y/t,z/t)

		bquarkmomentum.Boost(bCoM);

		long double cij, a, aprime;

		a = bquarkmomentum.Vect().Angle(motherLAB.Vect()); // get angle between v1 and v2
		// can it be computed from bCoM
		aprime = bquarkmomentum.Vect().Angle(-bCoM); // get angle between v1 and v2

		if (fabs(1 - aprime / a) < 0.01 || (a == aprime)) {
			cij = cos(a);

		} else {

			cout << "the cosine in the CoM does not match " << a << "!="
					<< aprime << endl;
			exit(1);
		}

		//	cout << "angolo e coseno fermione-W "<< a << " " << cij << endl;
		return cij;
	} else if (vectorstoadd.size() == 1) {

		cout
				<< "the cosine in the CoM can be computed only for two or more particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double cosij_cm(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of final state particles)

//cout<< vectorstoadd.size() <<endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 2) {

		long double cij, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; part[vect]->Status; vect++) {

			if (part[vect]->Status == 1) {
				sum = sum + vec_part[vect];
			}

			//		if (part[vectorstoadd[vect]]->PID == 11) {
			//			idvectminus = vect;
			//		}
			//		if (part[vectorstoadd[vect]]->PID == -13) {
			//			idvectplus = vect;
			//		}

		}

		Ycm = sum.Rapidity();

		//	Ycm = 0.5 * log((sum.E() + sum.Pz()) / (sum.E() - sum.Pz()));

		pplus = zBoost(vec_part[vectorstoadd[0]], Ycm);
		pminus = zBoost(vec_part[vectorstoadd[1]], Ycm);

		/*
		 pplus = vec_part[vectorstoadd[idvectminus]];
		 pminus = vec_part[vectorstoadd[idvectplus]];

		 Yplusprime = pplus.Rapidity() + Ycm;
		 Yminusprime = pminus.Rapidity() + Ycm;

		 pplus.SetE(pplus.Pt() * cosh(Yplusprime));
		 pplus.SetPz(pplus.Pt() * sinh(Yplusprime));
		 pminus.SetE(pminus.Pt() * cosh(Yminusprime));
		 pminus.SetPz(pminus.Pt() * sinh(Yminusprime));
		 */
		a = pplus.Vect().Angle(pminus.Vect()); // get angle between v1 and v2

		cij = cos(a);
		return cij;

	} else if (vectorstoadd.size() != 2) {

		cout << "the cosine in the CoM can be computed only for two particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double phi_gen(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of final state particles)

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

//	cout << vectorstoadd.size() << endl;

	if (vectorstoadd.size() == 4) {
		//		cout << "found 4 particles" << endl;

		long double cij, ctheta, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, kplus, kminus, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {

			if (part[vectorstoadd[vect]]->Status == -1) {

				sum = sum + vec_part[vectorstoadd[vect]];

				if (part[vectorstoadd[vect]]->PID < 0) {
					pminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					pplus = vec_part[vectorstoadd[vect]];
				}
			}

			if (part[vectorstoadd[vect]]->Status == 1) {
				if (part[vectorstoadd[vect]]->PID < 0) {
					kminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					kplus = vec_part[vectorstoadd[vect]];
				}
			}
		}

		/*  PHI RELATIVE */
		TVector3 eps1, eps2;
		TVector3 myaxis(1 / sqrt(3.), 1 / sqrt(3.), 1 / sqrt(3.));

		if ((pplus.Vect().Cross(pminus.Vect())).Mag()
				> 0.01 * (pplus.Vect().Mag() * pminus.Vect().Mag())) {

			eps1 = pplus.Vect().Cross(pminus.Vect());
			eps2 = sum.Vect().Cross(eps1);

		} else {

			eps1 = pplus.Vect().Cross(myaxis);
			eps2 = pplus.Vect().Cross(eps1);

		}

		float costheta, sintheta, cosphi, sinphi, phi;
		costheta = sum.Vect().Dot(kplus.Vect())
				/ (kplus.Vect().Mag() * sum.Vect().Mag());
		sintheta = sqrt(1 - costheta * costheta);

		cosphi = eps2.Dot(kplus.Vect())
				/ (sintheta * kplus.Vect().Mag() * eps2.Mag());
		sinphi = eps1.Dot(kplus.Vect())
				/ (sintheta * kplus.Vect().Mag() * eps1.Mag());

		phi = sign(sinphi) * acos(cosphi);

		return phi;

	} else if (vectorstoadd.size() != 4) {

		cout
				<< "the cos of the angle w.r.t incoming particle axis in the CoM can be computed only for four particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double cosij_cm2in(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of final state particles)

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

//	cout << vectorstoadd.size() << endl;

	if (vectorstoadd.size() == 4) {
		//		cout << "found 4 particles" << endl;

		long double cij, ctheta, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, kplus, kminus, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {

			if (part[vectorstoadd[vect]]->Status == -1) {

				sum = sum + vec_part[vectorstoadd[vect]];

				if (part[vectorstoadd[vect]]->PID < 0) {
					pminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					pplus = vec_part[vectorstoadd[vect]];
				}
			}

			if (part[vectorstoadd[vect]]->Status == 1) {
				if (part[vectorstoadd[vect]]->PID < 0) {
					kminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					kplus = vec_part[vectorstoadd[vect]];
				}
			}
		}

		/*  CoM polar angle */
		TVector3 b;
		b = sum.BoostVector();

		//		cout << "~~~~~~~~~~" << endl;
		//		cout << sum.Vect().Mag() << endl;

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);
		//		cout << temp.Vect().Mag() << endl;

		pplus.Boost(-b);
		pminus.Boost(-b);
		kplus.Boost(-b);
		kminus.Boost(-b);

		//		cout << kplus.Vect().Mag() << endl;
		//		cout << kminus.Vect().Mag() << endl;
		//		cout << (kplus + kminus).Vect().Mag() << endl;
		//		cout << (pplus + pminus).Vect().Mag() << endl;

		/*********************************/
		a = kplus.Vect().Angle(pplus.Vect()); // get angle between v1 and v2
		cij = cos(a);
		//		costheta = kplus.Vect().CosTheta();
		/*********************************/

		return cij;

	} else if (vectorstoadd.size() != 4) {

		cout
				<< "the cos of the angle w.r.t incoming particle axis in the CoM can be computed only for four particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double phi_cm2in(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of initial state particles)

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

//	cout << vectorstoadd.size() << endl;

	if (vectorstoadd.size() == 4) {
		//		cout << "found 4 particles" << endl;

		long double phi, ctheta, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, kplus, kminus, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {

			if (part[vectorstoadd[vect]]->Status == -1) {

				sum = sum + vec_part[vectorstoadd[vect]];

				if (part[vectorstoadd[vect]]->PID < 0) {
					pminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					pplus = vec_part[vectorstoadd[vect]];
				}
			}

			if (part[vectorstoadd[vect]]->Status == 1) {
				if (part[vectorstoadd[vect]]->PID < 0) {
					kminus = vec_part[vectorstoadd[vect]];
				}
				if (part[vectorstoadd[vect]]->PID > 0) {
					kplus = vec_part[vectorstoadd[vect]];
				}
			}
		}

		/*  CoM polar angle */
		TVector3 b;
		b = sum.BoostVector();
		//
		//		cout << "~~~~~~~~~~" << endl;
		//		cout << sum.Vect().Mag() << endl;

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);
		//		cout << temp.Vect().Mag() << endl;

		pplus.Boost(-b);
		pminus.Boost(-b);
		kplus.Boost(-b);
		kminus.Boost(-b);

		//		cout << kplus.Vect().Mag() << endl;
		//		cout << kminus.Vect().Mag() << endl;
		//		cout << (kplus + kminus).Vect().Mag() << endl;
		//		cout << (pplus + pminus).Vect().Mag() << endl;

		/*********************************/
		TVector3 eps1, eps2;

		eps1 = pplus.Vect().Cross(b);
		eps2 = pplus.Vect().Cross(eps1);

		/*********************************/

		float costheta, sintheta, cosphi, sinphi;
		costheta = pplus.Vect().Dot(kplus.Vect())
				/ (kplus.Vect().Mag() * pplus.Vect().Mag());
		sintheta = sqrt(1. - costheta * costheta);

		cosphi = eps2.Dot(kplus.Vect())
				/ (sintheta * kplus.Vect().Mag() * eps2.Mag());
		sinphi = eps1.Dot(kplus.Vect())
				/ (sintheta * kplus.Vect().Mag() * eps1.Mag());

		if (cosphi > 1. - 0.0001) {
			//			cout<<"pippo"<<endl;
			phi = 0;
		} else if (cosphi < -1. + 0.0001) {
			//			cout<<"pippo1"<<endl;
			phi = acos(-1.);
		} else {
			phi = sign(sinphi) * acos(cosphi);
		}
		return phi;

	} else if (vectorstoadd.size() != 4) {

		cout
				<< "the cos of the angle w.r.t incoming particle axis in the CoM can be computed only for four particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double costheta_cm2(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of final state particles)

//cout<< vectorstoadd.size() <<endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 2) {

		long double costheta, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {

			sum = sum + vec_part[vectorstoadd[vect]];

			//			hello();

			//		if (part[vectorstoadd[vect]]->PID == 11) {
			//			idvectminus = vect;
			//		}
			//		if (part[vectorstoadd[vect]]->PID == -13) {
			//			idvectplus = vect;
			//		}

		}

		TVector3 b;

		b = sum.BoostVector();
		//		cout << Ycm << endl;
		//	Ycm = 0.5 * log((sum.E() + sum.Pz()) / (sum.E() - sum.Pz()));

		//		cout << "~~~~~~~~~~" << endl;
		//		cout << sum.Vect().Mag() << endl;
		//		//			cout << (vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]]).Vect().Mag() << endl;
		//			cout << BoostZ(sum, Ycm).Vect().Mag() << endl;
		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);
		//		cout << temp.Vect().Mag() << endl;
		//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

		pplus = vec_part[vectorstoadd[0]];
		pplus.Boost(-b);
		pminus = vec_part[vectorstoadd[1]];
		pminus.Boost(-b);

		//		cout << pplus.Vect().Mag() << endl;
		//		cout << pminus.Vect().Mag() << endl;
		//		cout << (pplus + pminus).Vect().Mag() << endl;

		if (part[vectorstoadd[0]]->PID < 0) {
			costheta = pplus.Vect().CosTheta();
			if (isnan(costheta)) {
				costheta = vec_part[vectorstoadd[0]].Vect().CosTheta();

				cout << "nan found" << endl;
				cout << " now is " << costheta << endl;
			}
		}

		if (part[vectorstoadd[0]]->PID > 0) {
			costheta = pminus.Vect().CosTheta();
			if (isnan(costheta)) {
				costheta = vec_part[vectorstoadd[1]].Vect().CosTheta();
				cout << "nan found" << endl;
				cout << " now is " << costheta << endl;
			}

		}

		/*
		 pplus = vec_part[vectorstoadd[idvectminus]];
		 pminus = vec_part[vectorstoadd[idvectplus]];

		 Yplusprime = pplus.Rapidity() + Ycm;
		 Yminusprime = pminus.Rapidity() + Ycm;

		 pplus.SetE(pplus.Pt() * cosh(Yplusprime));
		 pplus.SetPz(pplus.Pt() * sinh(Yplusprime));
		 pminus.SetE(pminus.Pt() * cosh(Yminusprime));
		 pminus.SetPz(pminus.Pt() * sinh(Yminusprime));
		 */

		//		printlv(vec_part[vectorstoadd[0]]);
		//		printlv(vec_part[vectorstoadd[1]]);
		if (isnan(costheta)) {
			cout << "nan found" << endl;
			cout << " " << costheta << endl;
			printlv(vec_part[vectorstoadd[0]]);
			printlv(vec_part[vectorstoadd[1]]);
		}

		cout << " " << costheta << endl;

		return costheta;

	} else if (vectorstoadd.size() != 2) {

		cout
				<< "the cos of the angle w.r.t Z axis in the CoM can be computed only for two particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double costheta_cmZ2(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// compute the angle of the cosine between two particles in the center of mass frame
//(finding the frame summing all the momenta of final state particles)

//cout<< vectorstoadd.size() <<endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	if (vectorstoadd.size() == 2) {

		long double costheta, a, Ycm, Yplusprime, Yminusprime;

		TLorentzVector sum, pplus, pminus;
		int idvectplus, idvectminus;

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {

			sum = sum + vec_part[vectorstoadd[vect]];

			//			hello();

			//		if (part[vectorstoadd[vect]]->PID == 11) {
			//			idvectminus = vect;
			//		}
			//		if (part[vectorstoadd[vect]]->PID == -13) {
			//			idvectplus = vect;
			//		}

		}

		Ycm = sum.Rapidity();
		//		cout << Ycm << endl;
		//	Ycm = 0.5 * log((sum.E() + sum.Pz()) / (sum.E() - sum.Pz()));

		if (part[vectorstoadd[0]]->PID < 0) {
			//			cout << "~~~~~~~~~~" << endl;
			//			cout << sum.Vect().Mag() << endl;
			//			cout << (vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]]).Vect().Mag() << endl;
			//			cout << BoostZ(sum, Ycm).Vect().Mag() << endl;
			//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

			pplus = BoostZ(vec_part[vectorstoadd[0]], Ycm);
			pminus = BoostZ(vec_part[vectorstoadd[1]], Ycm);

			//			cout << pplus.Vect().Mag() << endl;
			//			cout << pminus.Vect().Mag() << endl;
			//			cout << (pplus + pminus).Vect().Mag() << endl;

			costheta = pplus.Vect().CosTheta();
			if (isnan(costheta)) {
				costheta = vec_part[vectorstoadd[0]].Vect().CosTheta();

				cout << "nan found" << endl;
				cout << " now is " << costheta << endl;
			}
		}

		if (part[vectorstoadd[0]]->PID > 0) {
			//			cout<< " PID neg"<<part[0]->PID <<endl;

			pplus = zBoost(vec_part[vectorstoadd[1]], Ycm);
			pminus = zBoost(vec_part[vectorstoadd[0]], Ycm);

			costheta = pplus.Vect().CosTheta();
			if (isnan(costheta)) {
				costheta = vec_part[vectorstoadd[1]].Vect().CosTheta();
				cout << "nan found" << endl;
				cout << " now is " << costheta << endl;
			}

		}

		/*
		 pplus = vec_part[vectorstoadd[idvectminus]];
		 pminus = vec_part[vectorstoadd[idvectplus]];

		 Yplusprime = pplus.Rapidity() + Ycm;
		 Yminusprime = pminus.Rapidity() + Ycm;

		 pplus.SetE(pplus.Pt() * cosh(Yplusprime));
		 pplus.SetPz(pplus.Pt() * sinh(Yplusprime));
		 pminus.SetE(pminus.Pt() * cosh(Yminusprime));
		 pminus.SetPz(pminus.Pt() * sinh(Yminusprime));
		 */

		//		printlv(vec_part[vectorstoadd[0]]);
		//		printlv(vec_part[vectorstoadd[1]]);
		if (isnan(costheta)) {
			cout << "nan found" << endl;
			cout << " " << costheta << endl;
			printlv(vec_part[vectorstoadd[0]]);
			printlv(vec_part[vectorstoadd[1]]);
		}

		//		cout << " " << costheta << endl;
		return costheta;

	} else if (vectorstoadd.size() != 2) {

		cout
				<< "the cos of the angle w.r.t Z axis in the CoM can be computed only for two particles"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double pTmaxOverMinvofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M, res;
	TLorentzVector sum, particle1;
	vector < TLorentzVector > fv;
	vector<float> pTs;

	int vect = 0;
//int n4vectors;
//n4vectors=sizeof(vec_part)/sizeof(vec_part[0]);
//cout<<"vettori letti"<<n4vectors<<endl;
//if ( size.vectorstoadd.size()==n4vectors ) {
	while (vect < vectorstoadd.size()) {

		sum = sum + vec_part[vectorstoadd[vect]];

		pTs.push_back(vec_part[vectorstoadd[vect]].Pt());
		fv.push_back(vec_part[vectorstoadd[vect]]);

		vect++;
	}

	M = sqrt(fabs(sum * sum));
//cout<<M<<endl;

	vector<int> idjets;
	idjets = MaxKElementsInAVector(pTs, 1);

//PrintVector(std::cout, idjets, " ");
//cout << endl;
//	PrintVector(std::cout, pTsJ, " ");
//	cout << endl;

	particle1 = fv[idjets[0]];

	res = particle1.Pt() / M;

//		cout << "The jet pT1 and pT2 are " << j1.Pt() << " " << j2.Pt() << endl;

	return res;
}

//****************************************************************************
extern "C" long double minvofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M;
	TLorentzVector sum;
	int vect = 0;
//int n4vectors;
//n4vectors=sizeof(vec_part)/sizeof(vec_part[0]);
//cout<<"vettori letti"<<n4vectors<<endl;
//if ( size.vectorstoadd.size()==n4vectors ) {
	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}
	M = sqrt(fabs(sum * sum));
//cout<<M<<endl;
	return M;
}

//****************************************************************************
extern "C" long double MpT(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M, result;
	TLorentzVector sum;
	int vect = 0;
//int n4vectors;
//n4vectors=sizeof(vec_part)/sizeof(vec_part[0]);
//cout<<"vettori letti"<<n4vectors<<endl;
//if ( size.vectorstoadd.size()==n4vectors ) {
	while (vect < vectorstoadd.size()) {
		sum = sum + vec_part[vectorstoadd[vect]];
		vect++;
	}
	M = sqrt(fabs(sum * sum));
//cout<<M<<endl;

	long double pTofn;

	vect = 0;
	TLorentzVector sump;

	while (vect < vectorstoadd.size()) {
		sump = sump + vec_part[vectorstoadd[vect]];
		vect++;
	}

	//printlv(sum);

	pTofn = sump.Pt();

	result = sqrt(M * M + pTofn * pTofn);
	return result;
}

//****************************************************************************
extern "C" long double err_pL_fwd_parton(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// err_px_fwd_parton 91 1
	long double result, diffp, diffm;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	TLorentzVector fwdp, fwdm, partm, partp;
	bool fwdp_exists, fwdm_exists, partonp_exists, partonm_exists;
	fwdp_exists = false;
	fwdm_exists = false;
	partonp_exists = false;
	partonm_exists = false;
	long double mindistancep, mindistancem;
	mindistancep = 1000;
	mindistancem = 1000;
	cout << "this event " << endl;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		int thispid;
		thispid = part[vectorstoadd[vect]]->PID;
		if (thispid == 91) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				fwdp = vec_part[vectorstoadd[vect]];
				fwdp_exists = true;
			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				fwdm = vec_part[vectorstoadd[vect]];
				fwdm_exists = true;
			}

		}
		// in an event with more hadron than just the forward jets
		// is necessary to find which parton level parton
		// has to be considered for the comparison
		//
		// i take the orginal parton that is closer in deltaR
		// deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(vec_part[vectorstoadd[1]].Vect());

		if (is_a_parton(thispid) == true) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				long double distancep;
				distancep = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdp.Vect()));
				cout << " this distance+ is " << distancep << endl;
				if (distancep < mindistancep) {
					cout << "new closest distance is " << distancep << endl;
					mindistancep = distancep;
					partp = vec_part[vectorstoadd[vect]];
					partonp_exists = true;
				}

			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				long double distancem;
				distancem = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdm.Vect()));
				if (distancem < mindistancem) {
					mindistancem = distancem;
					partm = vec_part[vectorstoadd[vect]];
					partonm_exists = true;
				}

			}

		}

	}

// now the loop on the particles is over
// and the momenta are assigned
// I have only to work out the result
	if (partonp_exists == true && fwdp_exists == true) {
		result = (fwdp.Pz() - partp.Pz()) / fwdp.Pz();
	} else {
		result = -10;
	}

//	if (partonm_exists == true && fwdm_exists == true) {
//		result = (fwdm.Px() - partm.Px()) / fwdm.Px();
//	} else {
//		result = -10;
//	}

	return result;
}

//****************************************************************************
extern "C" long double err_pT_fwd_parton(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// err_px_fwd_parton 91 1
	long double result, diffp, diffm;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	TLorentzVector fwdp, fwdm, partm, partp;
	bool fwdp_exists, fwdm_exists, partonp_exists, partonm_exists;
	fwdp_exists = false;
	fwdm_exists = false;
	partonp_exists = false;
	partonm_exists = false;
	long double mindistancep, mindistancem;
	mindistancep = 1000;
	mindistancem = 1000;
	cout << "this event " << endl;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		int thispid;
		thispid = part[vectorstoadd[vect]]->PID;
		if (thispid == 91) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				fwdp = vec_part[vectorstoadd[vect]];
				fwdp_exists = true;
			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				fwdm = vec_part[vectorstoadd[vect]];
				fwdm_exists = true;
			}

		}
		// in an event with more hadron than just the forward jets
		// is necessary to find which parton level parton
		// has to be considered for the comparison
		//
		// i take the orginal parton that is closer in deltaR
		// deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(vec_part[vectorstoadd[1]].Vect());

		if (is_a_parton(thispid) == true) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				long double distancep;
				distancep = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdp.Vect()));
				cout << " this distance+ is " << distancep << endl;
				if (distancep < mindistancep) {
					cout << "new closest distance is " << distancep << endl;
					mindistancep = distancep;
					partp = vec_part[vectorstoadd[vect]];
					partonp_exists = true;
				}

			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				long double distancem;
				distancem = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdm.Vect()));
				if (distancem < mindistancem) {
					mindistancem = distancem;
					partm = vec_part[vectorstoadd[vect]];
					partonm_exists = true;
				}

			}

		}

	}

// now the loop on the particles is ovesr
// and the momenta are assigned
// I have only to work out the result
	if (partonp_exists == true && fwdp_exists == true) {
		result = (fwdp.Pt() - partp.Pt()) / fwdp.Pt();
	} else {
		result = -10;
	}

//	if (partonm_exists == true && fwdm_exists == true) {
//		result = (fwdm.Px() - partm.Px()) / fwdm.Px();
//	} else {
//		result = -10;
//	}

	return result;
}

//****************************************************************************
extern "C" long double err_py_fwd_parton(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// err_px_fwd_parton 91 1
	long double result, diffp, diffm;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	TLorentzVector fwdp, fwdm, partm, partp;
	bool fwdp_exists, fwdm_exists, partonp_exists, partonm_exists;
	fwdp_exists = false;
	fwdm_exists = false;
	partonp_exists = false;
	partonm_exists = false;
	long double mindistancep, mindistancem;
	mindistancep = 1000;
	mindistancem = 1000;
	cout << "this event " << endl;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		int thispid;
		thispid = part[vectorstoadd[vect]]->PID;
		if (thispid == 91) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				fwdp = vec_part[vectorstoadd[vect]];
				fwdp_exists = true;
			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				fwdm = vec_part[vectorstoadd[vect]];
				fwdm_exists = true;
			}

		}
		// in an event with more hadron than just the forward jets
		// is necessary to find which parton level parton
		// has to be considered for the comparison
		//
		// i take the orginal parton that is closer in deltaR
		// deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(vec_part[vectorstoadd[1]].Vect());

		if (is_a_parton(thispid) == true) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				long double distancep;
				distancep = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdp.Vect()));
				cout << " this distance+ is " << distancep << endl;
				if (distancep < mindistancep) {
					cout << "new closest distance is " << distancep << endl;
					mindistancep = distancep;
					partp = vec_part[vectorstoadd[vect]];
					partonp_exists = true;
				}

			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				long double distancem;
				distancem = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdm.Vect()));
				if (distancem < mindistancem) {
					mindistancem = distancem;
					partm = vec_part[vectorstoadd[vect]];
					partonm_exists = true;
				}

			}

		}

	}

// now the loop on the particles is ovesr
// and the momenta are assigned
// I have only to work out the result
	if (partonp_exists == true && fwdp_exists == true) {
		result = (fwdp.Py() - partp.Py()) / fwdp.Py();
	} else {
		result = -10;
	}

//	if (partonm_exists == true && fwdm_exists == true) {
//		result = (fwdm.Px() - partm.Px()) / fwdm.Px();
//	} else {
//		result = -10;
//	}

	return result;
}

//****************************************************************************
extern "C" long double err_px_fwd_parton(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// err_px_fwd_parton 91 1
	long double result, diffp, diffm;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	TLorentzVector fwdp, fwdm, partm, partp;
	bool fwdp_exists, fwdm_exists, partonp_exists, partonm_exists;
	fwdp_exists = false;
	fwdm_exists = false;
	partonp_exists = false;
	partonm_exists = false;
	long double mindistancep, mindistancem;
	mindistancep = 1000;
	mindistancem = 1000;
	cout << "this event " << endl;
	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		int thispid;
		thispid = part[vectorstoadd[vect]]->PID;
		if (thispid == 91) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				fwdp = vec_part[vectorstoadd[vect]];
				fwdp_exists = true;
			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				fwdm = vec_part[vectorstoadd[vect]];
				fwdm_exists = true;
			}

		}
		// in an event with more hadron than just the forward jets
		// is necessary to find which parton level parton
		// has to be considered for the comparison
		//
		// i take the orginal parton that is closer in deltaR
		// deltaR = vec_part[vectorstoadd[0]].Vect().DeltaR(vec_part[vectorstoadd[1]].Vect());

		if (is_a_parton(thispid) == true) {
			if (vec_part[vectorstoadd[vect]].Pz() > 0) {
				long double distancep;
				distancep = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdp.Vect()));
				cout << " this distance+ is " << distancep << endl;
				if (distancep < mindistancep) {
					cout << "new closest distance is " << distancep << endl;
					mindistancep = distancep;
					partp = vec_part[vectorstoadd[vect]];
					partonp_exists = true;
				}

			}

			if (vec_part[vectorstoadd[vect]].Pz() < 0) {
				long double distancem;
				distancem = fabs(
						vec_part[vectorstoadd[vect]].Vect().DeltaPhi(
								fwdm.Vect()));
				if (distancem < mindistancem) {
					mindistancem = distancem;
					partm = vec_part[vectorstoadd[vect]];
					partonm_exists = true;
				}

			}

		}

	}

// now the loop on the particles is ovesr
// and the momenta are assigned
// I have only to work out the result
	if (partonp_exists == true && fwdp_exists == true) {
		result = (fwdp.Px() - partp.Px()) / fwdp.Px();
	} else {
		result = -10;
	}

//	if (partonm_exists == true && fwdm_exists == true) {
//		result = (fwdm.Px() - partm.Px()) / fwdm.Px();
//	} else {
//		result = -10;
//	}

	return result;
}

//****************************************************************************
extern "C" long double totalchargeofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double totcharge;
	totcharge = 0;

	for (int vect = 0; vect < vectorstoadd.size(); vect++) {

		int thispid;
		thispid = part[vectorstoadd[vect]]->PID;
		//cout<<thispid<<" "<<charge(thispid)<<endl;

		totcharge += charge(thispid);

	}

	cout << totcharge << endl;
	return totcharge;
}

//****************************************************************************
extern "C" long double bimass(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// this functions should only be used with PID=97 taken in pairs which return ordered pairs of fermins and antifermions

	long double M1, M2, res;
	TLorentzVector sum1, sum2;
	int vect = 0;

	if (matrixtoadd.size() == 2) {

		for (int vect = 0; vect < matrixtoadd[0].size(); vect++) {

			sum1 = sum1 + vec_part[matrixtoadd[0][vect]];

		}

		for (int vect = 0; vect < matrixtoadd[1].size(); vect++) {

			sum2 = sum2 + vec_part[matrixtoadd[1][vect]];

		}

		M1 = sqrt(fabs(sum1 * sum1));
		M2 = sqrt(fabs(sum2 * sum2));
		//		cout << M1 << " " << M2 << endl;
		//		return min(M2, M1) / max(M2, M1);

		//		cout << M1 << " " << M2 << endl;
		res = M1 + M2;
		cout << res << endl;
		return res;
	} else {

		cout << " bimass can be computed only on two jets at a time " << endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double PID(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	int pid;

	int vect = 0;

	if (vectorstoadd.size() == 1) {

		pid = part[vectorstoadd[0]]->PID;

		//cout<<M<<endl;
		return pid;
	} else {
		exit(1);
	}
}

//****************************************************************************
extern "C" long double mLHE(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M;

	int vect = 0;
//int n4vectors;
//n4vectors=sizeof(vec_part)/sizeof(vec_part[0]);
//cout<<"vettori letti"<<n4vectors<<endl;
//if ( size.vectorstoadd.size()==n4vectors ) {

	if (vectorstoadd.size() == 1) {

		M = part[vectorstoadd[0]]->M;

		//cout<<M<<endl;
		return M;
	} else {
		exit(1);
	}
}

//****************************************************************************
extern "C" long double HT1OverHT2_pair(int a, int b, int c, int d,
		vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
//****************************************************************************

	long double ht1, ht2;
	long double res;

//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	ht1 = vec_part[matrixtoadd[a][kp]].Pt() + vec_part[matrixtoadd[b][kp]].Pt();

	ht2 = vec_part[matrixtoadd[c][kp]].Pt() + vec_part[matrixtoadd[d][kp]].Pt();

	res = max(ht2, ht1) / min(ht1, ht2);

//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" long double deltaM_res_4into2(vector<int> ordering,
		vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
//****************************************************************************

	TLorentzVector pa, pb, pc, pd, sum1, sum2;
	long double deltares, M1, M2, detares, res;

//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	pa = vec_part[matrixtoadd[ordering[0]][kp]];
	pb = vec_part[matrixtoadd[ordering[1]][kp]];
	pc = vec_part[matrixtoadd[ordering[2]][kp]];
	pd = vec_part[matrixtoadd[ordering[3]][kp]];

	sum1 = pa + pb;
	sum2 = pc + pd;

	M1 = sqrt(fabs(sum1 * sum1));
	M2 = sqrt(fabs(sum2 * sum2));

	deltares = fabs((M1 - M2) / (M1 + M2));

	res = deltares;

//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
extern "C" long double deltaR_res_4into2(vector<int> ordering,
		vector<vector<int> > matrixtoadd, int kp,
		vector<TLorentzVector> vec_part) {
//****************************************************************************

	TLorentzVector pa, pb, pc, pd;
	long double deltares, dR1, dR2, detares, res;

//	cout<< matrixtoadd[0][0]<< " "<< matrixtoadd[1][0]<< " " << matrixtoadd[2][0]<<" "<< matrixtoadd[3][0]<<endl;

	pa = vec_part[matrixtoadd[ordering[0]][kp]];
	pb = vec_part[matrixtoadd[ordering[1]][kp]];
	pc = vec_part[matrixtoadd[ordering[2]][kp]];
	pd = vec_part[matrixtoadd[ordering[3]][kp]];

	dR1 = pa.Vect().DeltaR(pb.Vect());
	dR2 = pc.Vect().DeltaR(pd.Vect());

	deltares = fabs(dR1 - 1) + fabs(dR2 - 1);

	res = deltares;

//cout<< res.first << " "<< res.second << endl;

	return res;

}

//****************************************************************************
vector<struct CompositeParticle> FourIntoTwoByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	vector<long double> distance;
	distance.resize(0);

	TLorentzVector sum1, sum2, pa, pb, pc, pd;
//	vector < TLorentzVector > res;
//	pair < vector<TLorentzVector>, vector<TLorentzVector> > res;
	vector<struct CompositeParticle> res;

	struct CompositeParticle composite1, composite2;

	int element_min;

	int ord1[4] = { 0, 1, 2, 3 };
	int ord2[4] = { 0, 2, 1, 3 };
	int ord3[4] = { 0, 3, 2, 1 };

//int myints[] = { 16, 2, 77, 29 };
	vector<int> ordering1(ord1, ord1 + sizeof(ord1) / sizeof(int));
	vector<int> ordering2(ord2, ord2 + sizeof(ord2) / sizeof(int));
	vector<int> ordering3(ord3, ord3 + sizeof(ord3) / sizeof(int));

	vector < vector<int> > quartets;

	quartets.push_back(ordering1);
	quartets.push_back(ordering2);
	quartets.push_back(ordering3);

	if (matrixtoadd.size() == 4 && quartets.size() == 3) {

		for (int q = 0; q < quartets.size(); q++) {
			distance.push_back(
					deltaM_res_4into2(quartets[q], matrixtoadd, 0, vec_part)
//					deltaR_res_4into2(quartets[q], matrixtoadd, 0, vec_part)
							);

		}

		element_min = element_minimum_vector(distance);

		vector<int> bestquartet;
		bestquartet = quartets[element_min];

		pa = vec_part[matrixtoadd[bestquartet[0]][0]];
		pb = vec_part[matrixtoadd[bestquartet[1]][0]];
		pc = vec_part[matrixtoadd[bestquartet[2]][0]];
		pd = vec_part[matrixtoadd[bestquartet[3]][0]];

		vector<TLorentzVector> Vab, Vcd;

		Vab.push_back(pa);
		Vab.push_back(pb);

		Vcd.push_back(pc);
		Vcd.push_back(pd);

		sum1 = pa + pb;
		sum2 = pc + pd;

		composite1.Particle = sum1;
		composite1.Components = Vab;

		composite2.Particle = sum2;
		composite2.Components = Vcd;

		res.push_back(composite1);
		res.push_back(composite2);

		return res;

	} /* if 4 */else {
		cout << "not 4 particles in the reduction 4->2" << endl;
		exit(1);
	}

} //break in two resonances

//****************************************************************************
vector<struct CompositeParticle> FourIntoTwoByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	vector<long double> distance;
	distance.resize(0);

	TLorentzVector sum1, sum2, pa, pb, pc, pd;
//	vector < TLorentzVector > res;
//	pair < vector<TLorentzVector>, vector<TLorentzVector> > res;
	vector<struct CompositeParticle> res;

	struct CompositeParticle composite1, composite2;

	int element_min;

	int ord1[4] = { 0, 1, 2, 3 };
	int ord2[4] = { 0, 2, 1, 3 };
	int ord3[4] = { 0, 3, 2, 1 };

//int myints[] = { 16, 2, 77, 29 };
	vector<int> ordering1(ord1, ord1 + sizeof(ord1) / sizeof(int));
	vector<int> ordering2(ord2, ord2 + sizeof(ord2) / sizeof(int));
	vector<int> ordering3(ord3, ord3 + sizeof(ord3) / sizeof(int));

	vector < vector<int> > quartets;

	quartets.push_back(ordering1);
	quartets.push_back(ordering2);
	quartets.push_back(ordering3);

	if (matrixtoadd.size() == 4 && quartets.size() == 3) {

		for (int q = 0; q < quartets.size(); q++) {
			distance.push_back(
					deltaR_res_4into2(quartets[q], matrixtoadd, 0, vec_part));

		}

		element_min = element_minimum_vector(distance);

		vector<int> bestquartet;
		bestquartet = quartets[element_min];

		pa = vec_part[matrixtoadd[bestquartet[0]][0]];
		pb = vec_part[matrixtoadd[bestquartet[1]][0]];
		pc = vec_part[matrixtoadd[bestquartet[2]][0]];
		pd = vec_part[matrixtoadd[bestquartet[3]][0]];

		vector<TLorentzVector> Vab, Vcd;

		Vab.push_back(pa);
		Vab.push_back(pb);

		Vcd.push_back(pc);
		Vcd.push_back(pd);

		sum1 = pa + pb;
		sum2 = pc + pd;

		composite1.Particle = sum1;
		composite1.Components = Vab;

		composite2.Particle = sum2;
		composite2.Components = Vcd;

		res.push_back(composite1);
		res.push_back(composite2);

		return res;

	} /* if 4 */else {
		cout << "not 4 particles in the reduction 4->2" << endl;
		exit(1);
	}

} //break in two resonances

//****************************************************************************
extern "C" long double delta_eta_dijet_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, eta1, eta2, eta3, eta4;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		eta1 = ((ResonanceCandidates[0].Components)[0]).Eta();
		eta2 = ((ResonanceCandidates[0].Components)[1]).Eta();

		eta3 = ((ResonanceCandidates[1].Components)[0]).Eta();
		eta4 = ((ResonanceCandidates[1].Components)[1]).Eta();

		//res = min(fabs(eta1 - eta2), fabs(eta3 - eta4));
		res = 0.5 * (fabs(eta1 - eta2) + fabs(eta3 - eta4));

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double CosThetaStar_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, y1, y2, p0, p, costheta;

	TLorentzVector resonance1, resonance2, sum, resonance1star, resonance2star;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		sum = resonance1 + resonance2;

		//  y1 = resonance1.Rapidity();
		//	y2 = resonance2.Rapidity();

		TVector3 b;
		b = sum.BoostVector();

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);

		//cout << temp.Px() << " " << temp.Py() << " " << temp.Pz() << endl;

		//		cout << temp.Vect().Mag() << endl;
		//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

		resonance1star = resonance1;
		resonance1star.Boost(-b);

		resonance2star = resonance2;
		resonance2star.Boost(-b);

//		costheta = resonance1star.Vect().CosTheta();
		//	cout << costheta << endl;
		costheta = resonance1star.Pz() / (resonance1star.P());
		//	cout << costheta << endl;

		if (costheta + resonance2star.Vect().CosTheta() < 0.001) {

			res = costheta;

			return res;

		} else {
			cout << "boost to CoM did not give two back2back resonances"
					<< endl;

			cout << resonance1star.Vect().CosTheta() << " "
					<< resonance2star.Vect().CosTheta() << " " << endl;

			exit(1);
		}

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double AbsCosThetaStar_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, y1, y2, p0, p, costheta;

	TLorentzVector resonance1, resonance2, sum, resonance1star, resonance2star;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		sum = resonance1 + resonance2;

		//  y1 = resonance1.Rapidity();
		//	y2 = resonance2.Rapidity();

		TVector3 b;
		b = sum.BoostVector();

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);

		//cout << temp.Px() << " " << temp.Py() << " " << temp.Pz() << endl;

		//		cout << temp.Vect().Mag() << endl;
		//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

		resonance1star = resonance1;
		resonance1star.Boost(-b);

		resonance2star = resonance2;
		resonance2star.Boost(-b);

//		costheta = resonance1star.Vect().CosTheta();
		//	cout << costheta << endl;
		costheta = resonance1star.Pz() / (resonance1star.P());
		//	cout << costheta << endl;

		if (costheta + resonance2star.Vect().CosTheta() < 0.001) {

			res = fabs(costheta);

			return res;

		} else {
			cout << "boost to CoM did not give two back2back resonances"
					<< endl;

			cout << resonance1star.Vect().CosTheta() << " "
					<< resonance2star.Vect().CosTheta() << " " << endl;

			exit(1);
		}

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double AbsCosThetaStar_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, y1, y2, p0, p, costheta;

	TLorentzVector resonance1, resonance2, sum, resonance1star, resonance2star;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		sum = resonance1 + resonance2;

		//  y1 = resonance1.Rapidity();
		//	y2 = resonance2.Rapidity();

		TVector3 b;
		b = sum.BoostVector();

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);

		//cout << temp.Px() << " " << temp.Py() << " " << temp.Pz() << endl;

		//		cout << temp.Vect().Mag() << endl;
		//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

		resonance1star = resonance1;
		resonance1star.Boost(-b);

		resonance2star = resonance2;
		resonance2star.Boost(-b);

//		costheta = resonance1star.Vect().CosTheta();
		//	cout << costheta << endl;
		costheta = resonance1star.Pz() / (resonance1star.P());
		//	cout << costheta << endl;

		if (costheta + resonance2star.Vect().CosTheta() < 0.001) {

			res = fabs(costheta);

			return res;

		} else {
			cout << "boost to CoM did not give two back2back resonances"
					<< endl;

			cout << resonance1star.Vect().CosTheta() << " "
					<< resonance2star.Vect().CosTheta() << " " << endl;

			exit(1);
		}

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double CosThetaStar_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, y1, y2, p0, p, costheta;

	TLorentzVector resonance1, resonance2, sum, resonance1star, resonance2star;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		sum = resonance1 + resonance2;

		//  y1 = resonance1.Rapidity();
		//	y2 = resonance2.Rapidity();

		TVector3 b;
		b = sum.BoostVector();

		TLorentzVector temp;
		temp = sum;
		temp.Boost(-b);

		//cout << temp.Px() << " " << temp.Py() << " " << temp.Pz() << endl;

		//		cout << temp.Vect().Mag() << endl;
		//			cout << BoostZ(vec_part[vectorstoadd[0]] + vec_part[vectorstoadd[1]], Ycm).Vect().Mag() << endl;

		resonance1star = resonance1;
		resonance1star.Boost(-b);

		resonance2star = resonance2;
		resonance2star.Boost(-b);

//		costheta = resonance1star.Vect().CosTheta();
		//	cout << costheta << endl;
		costheta = resonance1star.Pz() / (resonance1star.P());
		//	cout << costheta << endl;

		if (costheta + resonance2star.Vect().CosTheta() < 0.001) {

			res = costheta;

			return res;

		} else {
			cout << "boost to CoM did not give two back2back resonances"
					<< endl;

			cout << resonance1star.Vect().CosTheta() << " "
					<< resonance2star.Vect().CosTheta() << " " << endl;

			exit(1);
		}

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_eta_max_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, eta1, eta2, eta3, eta4;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		eta1 = ((ResonanceCandidates[0].Components)[0]).Eta();
		eta2 = ((ResonanceCandidates[0].Components)[1]).Eta();

		eta3 = ((ResonanceCandidates[1].Components)[0]).Eta();
		eta4 = ((ResonanceCandidates[1].Components)[1]).Eta();

		//res = min(fabs(eta1 - eta2), fabs(eta3 - eta4));
		res = max(fabs(eta1 - eta2), fabs(eta3 - eta4));
		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_eta_min_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, eta1, eta2, eta3, eta4;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		eta1 = ((ResonanceCandidates[0].Components)[0]).Eta();
		eta2 = ((ResonanceCandidates[0].Components)[1]).Eta();

		eta3 = ((ResonanceCandidates[1].Components)[0]).Eta();
		eta4 = ((ResonanceCandidates[1].Components)[1]).Eta();

		//res = min(fabs(eta1 - eta2), fabs(eta3 - eta4));
		res = min(fabs(eta1 - eta2), fabs(eta3 - eta4));
		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_eta_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, eta1, eta2, eta3, eta4;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		eta1 = ((ResonanceCandidates[0].Components)[0]).Eta();
		eta2 = ((ResonanceCandidates[0].Components)[1]).Eta();

		eta3 = ((ResonanceCandidates[1].Components)[0]).Eta();
		eta4 = ((ResonanceCandidates[1].Components)[1]).Eta();

		//res = min(fabs(eta1 - eta2), fabs(eta3 - eta4));
		res = 0.5 * (fabs(eta1 - eta2) + fabs(eta3 - eta4));
		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_M_pseudoresonances_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, M1, M2;

	TLorentzVector resonance1, resonance2;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		M1 = sqrt(fabs(resonance1 * resonance1));
		M2 = sqrt(fabs(resonance2 * resonance2));

		res = fabs(M1 - M2) / (M1 + M2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double Mavg_pseudoresonances_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, M1, M2;

	TLorentzVector resonance1, resonance2;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		M1 = sqrt(fabs(resonance1 * resonance1));
		M2 = sqrt(fabs(resonance2 * resonance2));

		res = 0.5 * (M1 + M2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double Mavg_pseudoresonances_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, M1, M2;

	TLorentzVector resonance1, resonance2;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		M1 = sqrt(fabs(resonance1 * resonance1));
		M2 = sqrt(fabs(resonance2 * resonance2));

		res = 0.5 * (M1 + M2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double deltaR_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, dR1, dR2;

	TLorentzVector pa, pb, pc, pd;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		dR1 = pa.Vect().DeltaR(pb.Vect());
		dR2 = pc.Vect().DeltaR(pd.Vect());

		res = 0.5 * (dR1 + dR2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double deltaRmax_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, dR1, dR2;

	TLorentzVector pa, pb, pc, pd;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		dR1 = pa.Vect().DeltaR(pb.Vect());
		dR2 = pc.Vect().DeltaR(pd.Vect());

		res = max(dR1, dR2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double deltaRmin_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, dR1, dR2;

	TLorentzVector pa, pb, pc, pd;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		dR1 = pa.Vect().DeltaR(pb.Vect());
		dR2 = pc.Vect().DeltaR(pd.Vect());

		res = min(dR1, dR2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_phi_dijet_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, dphi1, dphi2;

	TLorentzVector pa, pb, pc, pd;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		dphi1 = fabs(pa.Vect().DeltaPhi(pb.Vect()));
		dphi2 = fabs(pc.Vect().DeltaPhi(pd.Vect()));

		res = 0.5 * (dphi1 + dphi2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double DeltaCMS_4Into2ByDeltaR(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, pta, ptb, ptc, ptd, mavg, dCMS1, dCMS2, dCMSmin, M1, M2;

	TLorentzVector pa, pb, pc, pd, resonance1, resonance2;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		M1 = sqrt(fabs(resonance1 * resonance1));
		M2 = sqrt(fabs(resonance2 * resonance2));

		mavg = (M1 + M2) / 2.;

		pta = pa.Pt();
		ptb = pb.Pt();
		ptc = pc.Pt();
		ptd = pd.Pt();

		dCMS1 = pta + ptb - mavg;
		dCMS2 = ptc + ptd - mavg;

		//dCMSavg = (dCMS1 + dCMS2) / 2.;
		dCMSmin = min(dCMS1, dCMS2);

		res = dCMSmin;

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_phi_dijet_4Into2ByDeltaR(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, dphi1, dphi2;

	TLorentzVector pa, pb, pc, pd;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaR(matrixtoadd, vec_part, part);

		pa = (ResonanceCandidates[0].Components)[0];
		pb = (ResonanceCandidates[0].Components)[1];
		pc = (ResonanceCandidates[1].Components)[0];
		pd = (ResonanceCandidates[1].Components)[1];

		dphi1 = fabs(pa.Vect().DeltaPhi(pb.Vect()));
		dphi2 = fabs(pc.Vect().DeltaPhi(pd.Vect()));

		res = 0.5 * (dphi1 + dphi2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double delta_M_pseudoresonances_4Into2ByDeltaM(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************

	long double res, M1, M2;

	TLorentzVector resonance1, resonance2;

	if (matrixtoadd.size() == 4) {

		vector<struct CompositeParticle> ResonanceCandidates;

		ResonanceCandidates = FourIntoTwoByDeltaM(matrixtoadd, vec_part, part);

		resonance1 = ResonanceCandidates[0].Particle;
		resonance2 = ResonanceCandidates[1].Particle;

		M1 = sqrt(fabs(resonance1 * resonance1));
		M2 = sqrt(fabs(resonance2 * resonance2));

		res = fabs(M1 - M2) / (M1 + M2);

		return res;

	} else {

		cout
				<< " any property of the best pairing can be computed only on 4 objects"
				<< endl;
		exit(1);

	}

}

//****************************************************************************
extern "C" long double HT_ratio_min(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printmatrix(matrixtoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	vector<long double> htratios;

	long double invratiomin, ratiomin;
	invratiomin = 0;
	ratiomin = 100.;

	if (matrixtoadd.size() == 4) {
		//HT12 = vec_part[vectorstoadd[0]].Pt()+ vec_part[vectorstoadd[1]].Pt();
		//HT34 = vec_part[vectorstoadd[2]].Pt()+ vec_part[vectorstoadd[3]].Pt();

		htratios.push_back(
				HT1OverHT2_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		htratios.push_back(
				HT1OverHT2_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		htratios.push_back(
				HT1OverHT2_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		for (int k = 0; k < htratios.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (1. / htratios[k] > invratiomin) {
				invratiomin = 1. / htratios[k];
			}

		}

		//cout << Mmin << endl;
		ratiomin = 1. / invratiomin;
		;

		ratiomin = min(htratios[0], min(htratios[1], htratios[2]));

		return ratiomin;

	} else {

		cout
				<< " the min{pairings}(H_T ratio) is coded for the case of 4 particles only"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double HT_ratio_max(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printmatrix(matrixtoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);
	vector<long double> htratios;

	long double ratiomax;
	ratiomax = 0;

	if (matrixtoadd.size() == 4) {
		//HT12 = vec_part[vectorstoadd[0]].Pt()+ vec_part[vectorstoadd[1]].Pt();
		//HT34 = vec_part[vectorstoadd[2]].Pt()+ vec_part[vectorstoadd[3]].Pt();

		htratios.push_back(
				HT1OverHT2_pair(0, 1, 2, 3, matrixtoadd, 0, vec_part));
		htratios.push_back(
				HT1OverHT2_pair(0, 2, 1, 3, matrixtoadd, 0, vec_part));
		htratios.push_back(
				HT1OverHT2_pair(0, 3, 2, 1, matrixtoadd, 0, vec_part));

		for (int k = 0; k < htratios.size(); k++) {
			//cout << delta_massa[k].first << " "<< delta_massa[k].second << endl;

			if (htratios[k] > ratiomax) {
				ratiomax = htratios[k];
			}

		}

		//cout << Mmin << endl;
		return ratiomax;

	} else {

		cout
				<< " the max{pairings}(H_T ratio) is coded for the case of 4 particles only"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double minvofnratio(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	long double M1, M2;
	TLorentzVector sum1, sum2;
	int vect = 0;
	if (matrixtoadd.size() == 2) {
		for (int vect = 0; vect < matrixtoadd[0].size(); vect++) {

			sum1 = sum1 + vec_part[matrixtoadd[0][vect]];

		}

		for (int vect = 0; vect < matrixtoadd[1].size(); vect++) {

			sum2 = sum2 + vec_part[matrixtoadd[1][vect]];

		}

		M1 = sqrt(fabs(sum1 * sum1));
		M2 = sqrt(fabs(sum2 * sum2));
		//cout<<M<<endl;
		return min(M2, M1) / max(M2, M1);
	} else {

		cout << " ratio of invariant masses can be computed only on two objects"
				<< endl;
		exit(1);

	}
}
//****************************************************************************
extern "C" long double minvofn2screenPID(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// this functions should only be used with status -23 from our custom Herwig

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M1, M2;
	TLorentzVector sum1, sum2;
	int vect = 0;

	if (vectorstoadd.size() == 2) {

		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			if (part[vectorstoadd[vect]]->PID > 0) {
				sum1 = vec_part[vectorstoadd[vect]];
			}

			if (part[vectorstoadd[vect]]->PID < 0) {
				sum2 = vec_part[vectorstoadd[vect]];
			}

		}

		M1 = sqrt(fabs(sum1 * sum1));
		M2 = sqrt(fabs(sum2 * sum2));
		cout << M1 << " " << M2 << endl;
		return min(M2, M1) / max(M2, M1);
	} else {

		cout << " " << endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double minvofn2screen(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

// this functions should only be used with PID=97 taken in pairs which return ordered pairs of fermins and antifermions

	long double M1, M2;
	TLorentzVector sum1, sum2;
	int vect = 0;

	if (matrixtoadd.size() == 2) {

		for (int vect = 0; vect < matrixtoadd[0].size(); vect++) {

			sum1 = sum1 + vec_part[matrixtoadd[0][vect]];

		}

		for (int vect = 0; vect < matrixtoadd[1].size(); vect++) {

			sum2 = sum2 + vec_part[matrixtoadd[1][vect]];

		}

		M1 = sqrt(fabs(sum1 * sum1));
		M2 = sqrt(fabs(sum2 * sum2));
		cout << M1 << " " << M2 << endl;
		return min(M2, M1) / max(M2, M1);
	} else {

		cout << " " << endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double minvofnratiofixed(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the total invariant mass of N particles
// particle involved are those in the event with IDs listed in  the vector vectorstoadd
//printvector(vectorstoadd);

	long double M1, M2;
	TLorentzVector sum1, sum2;
	int vect = 0;
	if (matrixtoadd.size() == 2) {
		for (int vect = 0; vect < matrixtoadd[0].size(); vect++) {

			sum1 = sum1 + vec_part[matrixtoadd[0][vect]];

		}

		for (int vect = 0; vect < matrixtoadd[1].size(); vect++) {

			sum2 = sum2 + vec_part[matrixtoadd[1][vect]];

		}

		M1 = sqrt(fabs(sum1 * sum1));
		M2 = sqrt(fabs(sum2 * sum2));
		//cout<<M<<endl;
		return M1 / M2;
	} else {

		cout << " ratio of invariant masses can be computed only on two objects"
				<< endl;
		exit(1);

	}
}

//****************************************************************************
extern "C" long double kTmeasureofn(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************

	if (matrixtoadd.size() == 2) {

		long double kTmeasure;
		TLorentzVector sum1;
		TLorentzVector sum2;

		//int n4vectors;
		//n4vectors=sizeof(vec_part)/sizeof(vec_part[0]);
		//cout<<"vettori letti"<<n4vectors<<endl;
		//if ( size.vectorstoadd.size()==n4vectors ) {
		for (int vect1 = 0; vect1 < matrixtoadd[0].size(); vect1++) {
			sum1 = sum1 + vec_part[matrixtoadd[0][vect1]];

		}

		for (int vect2 = 0; vect2 < matrixtoadd[1].size(); vect2++) {
			sum2 = sum2 + vec_part[matrixtoadd[1][vect2]];

		}

		long double E1, E2, EminQ;
		E1 = sum1.E();
		E2 = sum2.E();
		EminQ = min(E1 * E1, E2 * E2);

		long double a;
		a = sum1.Vect().Angle(sum2.Vect()); // get angle between v1 and v2

		kTmeasure = sqrt(2 * EminQ * (1 - cos(a)));

		return kTmeasure;
	} else {
		cout << "kT measure can be evaluated only between a pair of clusters"
				<< endl;
	}

}

//****************************************************************************
extern "C" long double t(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M;
	TLorentzVector sum, pf, pi;
	int vect = 0;
//printvector(vectorstoadd);

	if (vectorstoadd.size() == 2) {

		while (vect < vectorstoadd.size()) {
			if (part[vectorstoadd[vect]]->Status == 1) {
				pf = vec_part[vectorstoadd[vect]];
				//		cout<<"found final state"<<endl;
			}

			if (part[vectorstoadd[vect]]->Status == -1) {
				pi = vec_part[vectorstoadd[vect]];
				//				cout<<"found initial state"<<endl;

			}

			vect++;
		}

		sum = pf - pi;

		M = sqrt(fabs(sum * sum));

		//		vec_part[icount].SetPxPyPzE(part[icount]->Px,part[icount]->Py,part[icount]->Pz,part[icount]->E);
		//cout <<"im=" <<M << endl;
		return M;
	} else if (vectorstoadd.size() != 2) {

		cout << "t can be computed only for two particles" << endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double minvof2diff(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M;
	TLorentzVector sum;
	int vect = 0;
//printvector(vectorstoadd);

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		//cout << pid1 << " " << pid2 << endl;

		if (pid1 != pid2) {
			while (vect < vectorstoadd.size()) {

				sum = sum + vec_part[vectorstoadd[vect]];
				vect++;
			}
		}
		M = sqrt(sum * sum);
		//		vec_part[icount].SetPxPyPzE(part[icount]->Px,part[icount]->Py,part[icount]->Pz,part[icount]->E);
		//cout <<"im=" <<M << endl;
		return M;
	} else if (vectorstoadd.size() != 2) {

		cout
				<< "the invariant mass of two particles with different PID can be computed only for two particles"
				<< endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double minvof2fromZ(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double M;
	TLorentzVector sum;
	int vect = 0;
//printvector(vectorstoadd);

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		//cout << pid1 << " " << pid2 << endl;

		if (pid1 = -pid2) {
			while (vect < vectorstoadd.size()) {

				sum = sum + vec_part[vectorstoadd[vect]];
				vect++;
			}
		}
		M = sqrt(sum * sum);
		//		vec_part[icount].SetPxPyPzE(part[icount]->Px,part[icount]->Py,part[icount]->Pz,part[icount]->E);
		//cout <<"im=" <<M << endl;
		return M;
	} else if (vectorstoadd.size() != 2) {

		cout
				<< "the invariant mass of two particles with different PID can be computed only for two particles"
				<< endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double delta_phi_hjj(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vector<int> vectj;
	int vecthiggs;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double result, etah, etajone, etajtwo;
	TLorentzVector pjets, phiggs;
	vector < TLorentzVector > pj;

	if (vectorstoadd.size() == 3) {

		//	cout << "hey" << endl;
		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			if (part[vectorstoadd[vect]]->PID == 25) {

				phiggs = vec_part[vectorstoadd[vect]];
				vecthiggs = vectorstoadd[vect];

			}
		}

		for (int vect = 0; (vect < vectorstoadd.size()); vect++) {
			if (part[vectorstoadd[vect]]->PID != 25) {

				pj.push_back(vec_part[vectorstoadd[vect]]);
				vectj.push_back(vectorstoadd[vect]);
			}
		}

//		cout << "check session" << endl;

//		int pidh = (part[vecthiggs]->PID);
//		int pida = (part[vectj[0]]->PID);
//		int pidb = (part[vectj[1]]->PID);
//		cout << pidh << " " << pida << " " << pidb << endl;

		pjets = pj[0] + pj[1];
		result = fabs(phiggs.Vect().DeltaPhi(pjets.Vect())); // get angle between v1 and v2
		return result;
	} else if (vectorstoadd.size() != 3) {

		cout
				<< "the delta_phi_hjj variable can be computed only for 3 particles"
				<< endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double ZofZeppenfeld(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vector<int> vectj;
	int vecthiggs;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double result, etah, etajone, etajtwo;
	TLorentzVector pjone, pjtwo, phiggs;
	vector < TLorentzVector > pj;

	if (vectorstoadd.size() == 3) {

		//	cout << "hey" << endl;
		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			if (part[vectorstoadd[vect]]->PID == 25) {

				phiggs = vec_part[vectorstoadd[vect]];
				vecthiggs = vectorstoadd[vect];

			}
		}

		for (int vect = 0; (vect < vectorstoadd.size()); vect++) {
			if (part[vectorstoadd[vect]]->PID != 25) {

				pj.push_back(vec_part[vectorstoadd[vect]]);
				vectj.push_back(vectorstoadd[vect]);
			}
		}

//		cout << "check session" << endl;

//		int pidh = (part[vecthiggs]->PID);
//		int pida = (part[vectj[0]]->PID);
//		int pidb = (part[vectj[1]]->PID);
//		cout << pidh << " " << pida << " " << pidb << endl;

		etah = phiggs.Eta();
		etajone = pj[0].Eta();
		etajtwo = pj[1].Eta();
		result = etah - 0.5 * (etajone + etajtwo);
		return result;
	} else if (vectorstoadd.size() != 3) {

		cout << "the Z variable can be computed only for 3 particles" << endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double ZyyofZeppenfeld(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vector<int> vectj, vectgamma;
	int vecthiggs;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double result, etah, etajone, etajtwo;
	TLorentzVector pjone, pjtwo, pyone, pytwo;

	TLorentzVector phiggs;
	vector < TLorentzVector > pj;
	vector < TLorentzVector > pgamma;

	if (vectorstoadd.size() == 4) {

		//	cout << "hey" << endl;
		for (int vect = 0; vect < vectorstoadd.size(); vect++) {
			if (part[vectorstoadd[vect]]->PID == 22) {

				pgamma.push_back(vec_part[vectorstoadd[vect]]);
				vectgamma.push_back(vectorstoadd[vect]);

			}
		}

		for (int vect = 0; (vect < vectorstoadd.size()); vect++) {
			if (part[vectorstoadd[vect]]->PID != 22) {

				pj.push_back(vec_part[vectorstoadd[vect]]);
				vectj.push_back(vectorstoadd[vect]);
			}
		}

//		cout << "check session" << endl;

//		int pidh = (part[vecthiggs]->PID);
//		int pida = (part[vectj[0]]->PID);
//		int pidb = (part[vectj[1]]->PID);
//		cout << pidh << " " << pida << " " << pidb << endl;

		phiggs = pgamma[0] + pgamma[1];

		etah = phiggs.Eta();
		etajone = pj[0].Eta();
		etajtwo = pj[1].Eta();
		result = etah - 0.5 * (etajone + etajtwo);
		return result;
	} else if (vectorstoadd.size() != 4) {

		cout << "the Z variable can be computed only for 4 particles" << endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double pdOverpu(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of 2 particles with different PID
// this check is useful when one has a class made of 11,-11 and want to reconstruct a Z that
// must be reconstructed from a 11,-11 and not from a 11,11 pair.

//cout << "hello" << vectorstoadd.size()<< endl;

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double result;
	TLorentzVector pup, pdown;

	if (vectorstoadd.size() == 2) {

		int pid1 = part[vectorstoadd[0]]->PID;
		int pid2 = part[vectorstoadd[1]]->PID;
		cout << pid1 << " " << pid2 << endl;

		if (pid1 == 1 && pid2 == 2) {

			pdown = vec_part[vectorstoadd[0]];
			pup = vec_part[vectorstoadd[1]];

			result = sqrt(
					pdown.Px() * pdown.Px() + pdown.Py() * pdown.Py()
							+ pdown.Pz() * pdown.Pz())
					/ sqrt(
							pup.Px() * pup.Px() + pup.Py() * pup.Py()
									+ pup.Pz() * pup.Pz());

		} else if (pid1 == 2 && pid2 == 1) {

			pdown = vec_part[vectorstoadd[1]];
			pup = vec_part[vectorstoadd[0]];

			result = sqrt(
					pdown.Px() * pdown.Px() + pdown.Py() * pdown.Py()
							+ pdown.Pz() * pdown.Pz())
					/ sqrt(
							pup.Px() * pup.Px() + pup.Py() * pup.Py()
									+ pup.Pz() * pup.Pz());

		} else {
			result = -1.;
		}
		cout << result << endl;
		return result;
	} else if (vectorstoadd.size() != 2) {

		cout
				<< "the ratio of two particles momenta can be computed only for two particles"
				<< endl;
		return 0;

	}
}

//****************************************************************************
extern "C" long double ratio_res(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of all the pair of particles formed with the provided input
// if one pair has the same mass than any other pair of different particles the event
// contains a pair of resonances of equal mass.
//	cout << "**********" << endl;
//printvector(vectorstoadd);

// form the pairs

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	vector < vector<vector<int> > > allthepairs;
	findkplets(allthepairs, 2, vectorstoadd.size());

//cout << allthepairs.size() << endl;// one. becuase I am asking only the pairs
//cout << allthepairs[0].size() << endl; // the same as the number of pairs in the set
//printmatrix(allthepairs[0]);
// take two elements in the container of the pairs

	for (int pair1 = 0; pair1 < allthepairs[0].size(); pair1++) {
		for (int pair2 = 0; pair2 < allthepairs[0].size(); pair2++) {

			// cout << "selected lines " << pair1 << " " << pair2 << endl;

			// check particles are not used twice
			if (repeting_particle(allthepairs[0][pair1], allthepairs[0][pair2])
					== 0) {

				TLorentzVector v1, v2, v3, v4, othersum, sum;
				long double m1, m2;

				v1 = vec_part[vectorstoadd[allthepairs[0][pair1][0] - 1]];
				v2 = vec_part[vectorstoadd[allthepairs[0][pair1][1] - 1]];
				v3 = vec_part[vectorstoadd[allthepairs[0][pair2][0] - 1]];
				v4 = vec_part[vectorstoadd[allthepairs[0][pair2][1] - 1]];

				sum = v1 + v2;
				m1 = sum.Mag();

				othersum = v3 + v4;
				m2 = othersum.Mag();

				// cout << m2 << " " << m1 << endl;

				//				if ((m2 / m1 > 0.75) && (m2 / m1 < 1.25)) {
				if (fabs(m2 - m1) < 50) {
					//puts("trovato");
					//cout << m2 << " " << m1 << endl;
					return 1;
				}

			}
		}
	}
	return 0;
}
//****************************************************************************
extern "C" long double ratio_res_prop(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of all the pair of particles formed with the provided input
// if one pair has the same mass than any other pair of different particles the event
// contains a pair of resonances of equal mass.
//	cout << "**********" << endl;
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

// form the pairs
	vector < vector<vector<int> > > allthepairs;
	findkplets(allthepairs, 2, vectorstoadd.size());

//cout << allthepairs.size() << endl;// one. becuase I am asking only the pairs
//cout << allthepairs[0].size() << endl; // the same as the number of pairs in the set
//printmatrix(allthepairs[0]);
// take two elements in the container of the pairs

	for (int pair1 = 0; pair1 < allthepairs[0].size(); pair1++) {
		for (int pair2 = 0; pair2 < allthepairs[0].size(); pair2++) {

			// cout << "selected lines " << pair1 << " " << pair2 << endl;

			// check particles are not used twice
			if (repeting_particle(allthepairs[0][pair1], allthepairs[0][pair2])
					== 0) {

				TLorentzVector v1, v2, v3, v4, othersum, sum;
				long double m1, m2, ratio;

				v1 = vec_part[vectorstoadd[allthepairs[0][pair1][0] - 1]];
				v2 = vec_part[vectorstoadd[allthepairs[0][pair1][1] - 1]];
				v3 = vec_part[vectorstoadd[allthepairs[0][pair2][0] - 1]];
				v4 = vec_part[vectorstoadd[allthepairs[0][pair2][1] - 1]];

				sum = v1 + v2;
				m1 = sum.Mag();

				othersum = v3 + v4;
				m2 = othersum.Mag();

				ratio = m2 / m1;
				if ((ratio > 0.75) && (ratio < 1.25)) {
					//				if ( fabs(m2 - m1) < 50 ) {
					//puts("trovato");
					//cout << m2 << " " << m1 << endl;
					return 1;
				}

			}
		}
	}
	return 0;
}

//****************************************************************************
extern "C" long double best_pair_res_prop(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of all the pair of particles formed with the provided input
// if one pair has the same mass than any other pair of different particles the event
// contains a pair of resonances of equal mass.
//	cout << "**********" << endl;
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

// form the pairs
	vector < vector<vector<int> > > allthepairs;
	findkplets(allthepairs, 2, vectorstoadd.size());

//cout << allthepairs.size() << endl;// one. becuase I am asking only the pairs
//cout << allthepairs[0].size() << endl; // the same as the number of pairs in the set
//printmatrix(allthepairs[0]);
// take two elements in the container of the pairs
	long double res;
	res = 0;
	float window = 0.25;

	for (int pair1 = 0; pair1 < allthepairs[0].size(); pair1++) {
		for (int pair2 = 0; pair2 < allthepairs[0].size(); pair2++) {

			// cout << "selected lines " << pair1 << " " << pair2 << endl;

			// check particles are not used twice
			if (repeting_particle(allthepairs[0][pair1], allthepairs[0][pair2])
					== 0) {

				TLorentzVector v1, v2, v3, v4, othersum, sum;
				long double m1, m2, ratio, bestratio;
				bestratio = 1 - window * 1.01;
				//				cout<< bestratio<<endl;
				v1 = vec_part[vectorstoadd[allthepairs[0][pair1][0] - 1]];
				v2 = vec_part[vectorstoadd[allthepairs[0][pair1][1] - 1]];
				v3 = vec_part[vectorstoadd[allthepairs[0][pair2][0] - 1]];
				v4 = vec_part[vectorstoadd[allthepairs[0][pair2][1] - 1]];

				sum = v1 + v2;
				m1 = sum.Mag();

				othersum = v3 + v4;
				m2 = othersum.Mag();

				ratio = m2 / m1;

				if ((ratio > 1 - window) && (ratio < 1 + window)) {
					//					cout<<ratio<<endl;
					if (fabs(bestratio - 1) > fabs(ratio - 1)) {
						bestratio = ratio;
						res = 0.5 * (m2 + m1);
					}

				} //if within the window

			} //if different particles
		}
	}
	return res; //if no good pairing has been found the result will be zero
}

//****************************************************************************
extern "C" long double ratio_resT_prop(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the  invariant mass of all the pair of particles formed with the provided input
// if one pair has the same mass than any other pair of different particles the event
// contains a pair of resonances of equal mass.
//	cout << "**********" << endl;
//printvector(vectorstoadd);

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

// form the pairs
	vector < vector<vector<int> > > allthepairs;
	findkplets(allthepairs, 2, vectorstoadd.size());

//cout << allthepairs.size() << endl;// one. becuase I am asking only the pairs
//cout << allthepairs[0].size() << endl; // the same as the number of pairs in the set
//printmatrix(allthepairs[0]);
// take two elements in the container of the pairs

	for (int pair1 = 0; pair1 < allthepairs[0].size(); pair1++) {
		for (int pair2 = 0; pair2 < allthepairs[0].size(); pair2++) {

			// cout << "selected lines " << pair1 << " " << pair2 << endl;

			// check particles are not used twice
			if (repeting_particle(allthepairs[0][pair1], allthepairs[0][pair2])
					== 0) {

				TLorentzVector v1, v2, v3, v4;
				long double mT1, mT2;

				v1 = vec_part[vectorstoadd[allthepairs[0][pair1][0] - 1]];
				v2 = vec_part[vectorstoadd[allthepairs[0][pair1][1] - 1]];
				v3 = vec_part[vectorstoadd[allthepairs[0][pair2][0] - 1]];
				v4 = vec_part[vectorstoadd[allthepairs[0][pair2][1] - 1]];

				long double cTij1, cTij2, deltaphi1, deltaphi2, E1i, E1j, E2i,
						E2j;

				deltaphi1 = v1.Vect().DeltaPhi(v2.Vect());

				cTij1 = cos(deltaphi1);
				E1i = v1.E();
				E1j = v2.E();
				mT1 = sqrt(2 * E1i * E1j * (1 - cTij1));

				deltaphi2 = v3.Vect().DeltaPhi(v4.Vect());

				cTij2 = cos(deltaphi2);
				E2i = v3.E();
				E2j = v4.E();
				mT2 = sqrt(2 * E2i * E2j * (1 - cTij2));

				if ((mT2 / mT1 > 0.9) && (mT2 / mT1 < 1.1)) {
					//				if ( fabs(m2 - m1) < 50 ) {
					//puts("trovato");
					//cout << m2 << " " << m1 << endl;
					return 1;
				}

			}
		}
	}
	return 0;
}

extern "C" long double transverse_sphericity(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the angle of the cosine brtween two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double mXX, mYY, mXY, ST, thisPx, thisPy, minEigen, maxEigen;

	mXX = 0;
	mYY = 0;
	mXY = 0;

	for (int part = 0; part < vectorstoadd.size(); part++) {
		thisPx = vec_part[vectorstoadd[part]].Px();
		thisPy = vec_part[vectorstoadd[part]].Py();
		mXX = mXX + thisPx * thisPx;
		mYY = mYY + thisPy * thisPy;
		mXY = mXY + thisPx * thisPy;
		//mYX=mYX+thisPx*thisPy;
	}

	minEigen = (mXX + mYY
			- sqrt(mXX * mXX + 4 * mXY * mXY - 2 * mXX * mYY + mYY * mYY)) / 2.;
	maxEigen = (mXX + mYY
			+ sqrt(mXX * mXX + 4 * mXY * mXY - 2 * mXX * mYY + mYY * mYY)) / 2.;

	ST = 2 * minEigen / (minEigen + maxEigen);

	return ST; //} else {

//}
}

extern "C" long double transverse_sphericity_ROOT(
		vector<vector<int> > matrixtoadd, vector<TLorentzVector> vec_part,
		vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the angle of the cosine brtween two particles in the lab frame

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double mXX, mYY, mXY, ST, thisPx, thisPy, minEigen, maxEigen;

	mXX = 0;
	mYY = 0;
	mXY = 0;

	for (int part = 0; part < vectorstoadd.size(); part++) {
		thisPx = vec_part[vectorstoadd[part]].Px();
		thisPy = vec_part[vectorstoadd[part]].Py();
		mXX = mXX + thisPx * thisPx;
		mYY = mYY + thisPy * thisPy;
		mXY = mXY + thisPx * thisPy;
		//mYX=mYX+thisPx*thisPy;
	}

	const int N = 2;

	double e[N * N] = {
		mXX, mXY,
		mXY, mYY
	};

TMatrixDSym 	m(N, e);
	TMatrixDSymEigen me(m);

	TVectorD eigenval = me.GetEigenValues();

// m.Print();

// cout << mXX << " " << mXY << endl;
// cout << mXY << " " << mYY << endl;

//minEigen = (mXX + mYY
//		- sqrt(mXX * mXX + 4 * mXY * mXY - 2 * mXX * mYY + mYY * mYY)) / 2.;
//maxEigen = (mXX + mYY
//	+ sqrt(mXX * mXX + 4 * mXY * mXY - 2 * mXX * mYY + mYY * mYY)) / 2.;

	minEigen = eigenval[1];
	maxEigen = eigenval[0];

//	eigenval.Print();
//cout << minEigen << " " << maxEigen << endl;
//cout << eigenval[1] << " " << eigenval[0] << endl;

	ST = 2 * minEigen / (minEigen + maxEigen);

	return ST; //} else {

//}
}

extern "C" long double sphericity(vector<vector<int> > matrixtoadd,
		vector<TLorentzVector> vec_part, vector<TRootLHEFParticle*> part) {
//****************************************************************************
// compute the sphericity as in J.D. Bjorken and S.J. Brodsky, Phys. Rev. D1 (1970) 1416

	vector<int> vectorstoadd;
	vectorstoadd = matrix2vector(matrixtoadd);

	long double mXX, mYY, mXY, mXZ, mYZ, mZZ, S3D, thisPx, thisPy, thisPz,
			minEigen, nexttominEigen, sumAbsP;

	mXX = 0;
	mYY = 0;
	mXY = 0;
	mXZ = 0;
	mYZ = 0;
	mZZ = 0;
	sumAbsP = 0;

	for (int part = 0; part < vectorstoadd.size(); part++) {
		thisPx = vec_part[vectorstoadd[part]].Px();
		thisPy = vec_part[vectorstoadd[part]].Py();
		thisPz = vec_part[vectorstoadd[part]].Pz();

		sumAbsP = sumAbsP + thisPx * thisPx + thisPy * thisPy + thisPz * thisPz;

		mXX = mXX + thisPx * thisPx;
		mYY = mYY + thisPy * thisPy;
		mXY = mXY + thisPx * thisPy;
		mXZ = mXZ + thisPx * thisPz;
		mYZ = mYZ + thisPy * thisPz;
		mZZ = mZZ + thisPz * thisPz;
		//mYX=mYX+thisPx*thisPy;
	}

	mXX = mXX / sumAbsP;
	mXY = mXY / sumAbsP;
	mYY = mYY / sumAbsP;

	mXZ = mXZ / sumAbsP;
	mYZ = mYZ / sumAbsP;
	mZZ = mZZ / sumAbsP;

	const int N = 3;

	double e[N * N] = {
		mXX, mXY, mXZ,
		mXY, mYY, mYZ,
		mXZ, mYZ, mZZ
	};

TMatrixDSym 	m(N, e);
	TMatrixDSymEigen me(m);

	TVectorD eigenval = me.GetEigenValues();

//m.Print();
//	eigenval.Print();

	minEigen = eigenval[2];
	nexttominEigen = eigenval[1];

	S3D = 3 / 2. * (minEigen + nexttominEigen);

	return S3D; //} else {

//}
}

#endif
